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Abstract 

In  this  paper  concepts  from  continuum  mechanics  are  used  to  define  geodesic  paths 
in  the  space  of  shapes,  where  shapes  are  implicitly  described  as  boundary  contours  of 
objects.  The  proposed  shape  metric  is  derived  from  a  continuum  mechanical  notion  of 
viscous  dissipation.  A  geodesic  path  is  defined  as  the  family  of  shapes  such  that  the 
total  amount  of  viscous  dissipation  caused  by  an  optimal  material  transport  along  the 
path  is  minimized.  The  approach  can  easily  be  generalized  to  shapes  given  as  segment 
contours  of  multi-labeled  images  and  to  geodesic  paths  between  partially  occluded  ob¬ 
jects.  The  proposed  computational  framework  for  finding  such  a  minimizer  is  based  on 
the  time  discretization  of  a  geodesic  path  as  a  sequence  of  pairwise  matching  problems, 
which  is  strictly  invariant  with  respect  to  rigid  body  motions  and  ensures  a  1-1  corre¬ 
spondence  along  the  induced  flow  in  shape  space.  When  decreasing  the  time  step  size, 
the  proposed  model  leads  to  the  minimization  of  the  actual  geodesic  length,  where  the 
Hessian  of  the  pairwise  matching  energy  reflects  the  chosen  Riemannian  metric  on  the 
underlying  shape  space.  If  the  constraint  of  pairwise  shape  correspondence  is  replaced 
by  the  volume  of  the  shape  mismatch  as  a  penalty  functional,  one  obtains  for  decreas¬ 
ing  time  step  size  an  optical  flow  term  controlling  the  transport  of  the  shape  by  the 
underlying  motion  field.  The  method  is  implemented  via  a  level  set  representation  of 
shapes,  and  a  finite  element  approximation  is  employed  as  spatial  discretization  both 
for  the  pairwise  matching  deformations  and  for  the  level  set  representations.  The  nu¬ 
merical  relaxation  of  the  energy  is  performed  via  an  efficient  multi-scale  procedure  in 
space  and  time.  Various  examples  for  2D  and  3D  shapes  underline  the  effectiveness 
and  robustness  of  the  proposed  approach. 


1  Introduction 

In  this  paper  we  investigate  the  close  link  between  abstract  geometry  on  the  infinite-dimen¬ 
sional  space  of  shapes  and  the  continuum  mechanical  view  of  shapes  as  boundary  contours 
of  physical  objects  in  order  to  define  geodesic  paths  and  distances  between  shapes  in  2D  and 
3D.  The  computation  of  shape  distances  and  geodesics  is  fundamental  for  problems  ranging 
from  computational  anatomy  to  object  recognition,  warping,  and  matching.  The  aim  is  to 
reliably  and  effectively  evaluate  distances  between  non-parametrized  geometric  shapes  of 
possibly  different  topology.  In  particular,  we  allow  shapes  to  consist  of  boundary  contours 
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Figure  1:  Time-discrete  geodesic  between  the  letters  A  and  B.  The  geodesic  distance  is 
measured  on  the  basis  of  viscous  dissipation  inside  the  objects  (color-coded  in  the  top  row 
from  blue,  low  dissipation,  to  red,  high  dissipation),  which  is  approximated  as  a  deformation 
energy  of  pairwise  1-1  deformations  between  consecutive  shapes  along  the  discrete  geodesic 
path.  Shapes  are  represented  via  level  set  functions,  whose  level  lines  are  texture-coded  in 
the  bottom  row. 

of  multiple  components  of  volumetric  objects.  The  underlying  Riemannian  metric  on  shape 
space  is  identified  with  physical  dissipation  (cf.  Fig.  1) — the  rate  at  which  mechanical  energy 
is  converted  into  heat  in  a  viscous  fluid  due  to  friction — accumulated  along  an  optimal 
transport  of  the  volumetric  objects  (cf.  [47]). 

We  simultaneously  address  the  following  major  challenges:  A  physically  sound  modeling 
of  the  geodesic  flow  of  shapes  given  as  boundary  contours  of  possibly  multi-component 
objects  on  a  void  background,  the  need  for  a  coarse  time  discretization  of  the  continuous 
geodesic  path,  and  a  numerically  effective  relaxation  of  the  resulting  time-  and  space-discrete 
variational  problem.  Addressing  these  challenges  leads  to  a  novel  formulation  for  discrete 
geodesic  paths  in  shape  space  that  is  based  on  solid  mathematical,  computational,  and 
physical  arguments  and  motivations. 

Different  from  the  pioneering  diffeomorphism  approach  by  Miller  et  al.  [35]  the  motion  held 
V  governing  the  flow  in  shape  space  vanishes  on  the  object  background,  and  the  accumulated 
physical  dissipation  is  a  quadratic  functional  depending  only  on  the  first  order  local  variation 
of  a  flow  field.  In  fact,  as  we  will  explain  in  a  separate  section  on  the  physical  background, 
the  dissipation  depends  only  on  the  symmetric  part  e[n]  =  ^(Dn^  +  T>v)  of  the  Jacobian  Vv 
of  the  motion  field  v,  and  under  the  additional  assumption  of  isotropy,  a  typical  model  for 
the  dissipation  is  given  by  Diss[n]  =  /g  diss[n]  dxdt  with  the  local  rate  of  dissipation 

diss[n]  =  —  (tre[n])^  +  /itr(e[n]^)  (1) 
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(cf.  [21]),  where  0(t)  describes  the  deformed  object.  The  outer  integral  accumulates  the 
dissipation  in  time  during  the  deformation  of  0(0)  into  0(1).  The  physical  variable  t  geo¬ 
metrically  represents  the  coordinate  along  the  path  in  shape  space. 

A  straightforward  time  discretization  of  a  geodesic  flow  would  neither  guarantee  local  rigid 
body  motion  invariance  for  the  time-discrete  problem  nor  a  1-1  mapping  between  objects 
at  consecutive  time  steps.  For  this  reason  we  present  a  time  discretization  which  is  based 
on  a  pairwise  matching  of  intermediate  shapes  that  correspond  to  subsequent  time  steps. 
In  fact,  such  a  discretization  of  a  path  as  concatenation  of  short  connecting  line  segments 
in  shape  space  between  consecutive  shapes  is  natural  with  regard  to  the  variational  defini¬ 
tion  of  a  geodesic.  It  also  underlies  for  instance  the  algorithm  by  Schmidt  et  al.  [37]  and 
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Figure  2:  Discrete  geodesics  between  a  straight  and  a  rolled  up  bar,  from  first  row  to  fourth 
row  based  on  1,  2,  4,  and  8  time  steps.  The  light  gray  shapes  in  the  first,  second,  and  third 
row  show  a  linear  interpolation  of  the  deformations  connecting  the  dark  gray  shapes.  The 
shapes  from  the  finest  time  discretization  are  overlayed  over  the  others  as  thin  black  lines. 
In  the  last  row  the  rate  of  viscous  dissipation  is  rendered  on  the  shape  domains  ...  ,0^ 
from  the  previous  row,  color-coded  as  Hi  ■. 


can  be  regarded  as  the  infinite-dimensional  counterpart  of  the  following  time  discretization 
for  a  geodesic  between  two  points  and  on  a  finite-dimensional  Riemannian  manifold: 
Consider  a  sequence  of  points  =  Sq,  Si,  . . . ,  sx  =  sb  connecting  two  fixed  points  sa 
and  Sb  and  minimize  dist^(sfc_i,  Sfc),  where  dist(-,  •)  is  a  suitable  approximation  of  the 

Riemannian  distance.  In  our  case  of  the  infinite-dimensional  shape  space,  dist^(-,  •)  will  be 
approximated  by  a  suitable  energy  of  the  matching  deformation  between  subsequent  shapes. 
In  particular,  we  will  employ  a  deformation  energy  from  the  class  of  so-called  polyconvex 
energies  [14]  to  ensure  both  exact  frame  indifference  (observer  independence  and  thus  rigid 
body  motion  invariance)  and  a  global  1-1  property.  Both  the  built-in  exact  frame  indiffer¬ 
ence  and  the  1-1  mapping  property  ensure  that  fairly  coarse  time  discretizations  already 
lead  to  an  accurate  approximation  of  geodesic  paths  (cf.  Fig.  2).  The  approach  is  inspired 
both  by  work  in  mechanics  [46]  and  in  geometry  [29] .  We  will  also  discuss  the  corresponding 
continuous  problem  when  the  time  discretization  step  vanishes. 

Careful  consideration  is  required  with  respect  to  the  effective  multi-scale  minimization  of 
the  time  discrete  path  length.  Already  in  the  case  of  low-dimensional  Riemannian  manifolds 
the  need  for  an  efficient  cascadic  coarse  to  fine  minimization  strategy  is  apparent.  To  give  a 
conceptual  sketch  of  the  proposed  algorithm  on  the  actual  shape  space,  Fig.  3  demonstrates 
the  proposed  procedure  in  the  case  of  considered  as  the  stereographic  projection  of  the 
two-dimensional  sphere,  which  already  illustrates  the  advantage  of  our  proposed  optimiza¬ 
tion  framework. 

The  organization  of  the  paper  is  as  follows.  Sections  1.1  and  1.2  respectively  give  a  brief 
introduction  to  the  continuum  mechanical  background  of  dissipation  in  viscous  fluid  trans- 
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Figure  3:  Different  refinement  levels  of  a  discrete  geodesic  [K  =  1,2,4,...,  256)  from  Johan¬ 
nesburg  to  New  York  in  the  stereographic  projection  (right)  and  backprojected  on  the  globe 
(left).  The  discrete  geodesic  for  a  given  K  minimizes  dist^(sfc_i,  s^),  where  the  are 

points  on  the  globe  (represented  by  the  black  dots  in  the  stereographic  projection)  and  Sq  and 
sk  correspond  to  Johannesburg  and  New  York,  respectively.  dist(sfc_i,  s^)  is  approximated 
by  measuring  the  length  of  the  segment  (sfc_i,Sfc)  in  the  stereographic  projection,  using  the 
stereographic  metric  at  the  segment  midpoint.  The  red  line  shows  the  discrete  geodesic  on 
the  finest  level.  A  single-level  nonlinear  Gauss-Seidel  relaxation  of  the  corresponding  energy 
on  the  finest  resolution  with  successive  relaxation  of  the  different  vertices  requires  over  10® 
elementary  relaxation  steps,  whereas  in  a  cascadic  energy  relaxation  scheme,  which  proceeds 
from  coarse  to  fine  resolution,  only  2579  of  these  elementary  minimization  steps  are  needed. 


port  and  discuss  related  work  on  shape  distances  and  geodesics  in  shape  space,  examining 
the  relation  to  physics.  Section  1.3  lists  the  key  contributions  of  our  approach.  Section  2  is 
devoted  to  the  proposed  variational  approach.  We  first  introduce  the  notion  of  time-discrete 
geodesics  in  Section  2.1,  prove  existence  under  suitable  assumptions  in  in  Section  2.2,  and 
we  present  a  relaxed  formulation  in  Section  2.3.  Then,  in  Section  2.4  we  present  the  actual 
viscous  fluid  model  for  geodesics  in  shape  space  and  establish  it  as  the  limit  model  of  our  time 
discretization  for  vanishing  time  step  size  in  Section  2.5.  Section  3  introduces  the  correspond¬ 
ing  numerical  algorithm,  wich  is  based  on  a  regularized  level  set  approximation  as  described 
in  Section  3.1  and  the  space  discretization  via  finite  elements  as  detailed  in  Section  3.2.  A 
sketch  of  the  proposed  overall  multi-scale  algorithm  is  provided  in  Section  3.3.  Section  4  is 
devoted  to  the  computational  results  and  various  applications,  including  geodesics  in  2D  and 
3D,  shapes  as  boundary  contours  of  multi-labeled  objects,  applications  to  shape  statistics, 
and  an  illustrative  analysis  of  parts  of  the  global  shape  space  structure.  Finally,  in  Section  5 
we  draw  conclusions  and  describe  prospective  research  directions. 

1.1  The  physical  background  revisited 

Our  approach  relies  on  a  close  link  between  geodesics  in  shape  space  and  the  continuum 
mechanics  of  viscous  fluid  transport.  Therefore,  we  will  here  review  the  fundamental  concept 
of  viscous  dissipation  in  a  Newtonian  fluid.  The  section  is  intended  for  readers  less  familiar 
with  this  topic  and  can  be  skipped  otherwise. 

Even  though  fluids  are  composed  of  molecules,  based  on  the  common  continuum  as¬ 
sumption  one  studies  the  macroscopic  behavior  of  a  fluid  via  governing  partial  differential 
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Figure  4:  A  linear  velocity  profile  produces  a  pure  horizontal  shear  stress. 


equations  which  describe  the  transport  of  fluid  material.  Here,  viscosity  describes  the  internal 
resistance  in  a  fluid  and  may  be  thought  of  as  a  macroscopic  measure  of  the  friction  between 
fluid  particles.  As  an  example,  the  viscosity  of  honey  is  significantly  larger  than  that  of 
water.  Mathematically,  the  friction  is  described  in  terms  of  the  stress  tensor  a  = 
whose  entries  describe  a  force  per  area  element.  By  definition,  cTjj  is  the  force  component 
along  the  Ah  coordinate  direction  acting  on  the  area  element  with  a  normal  pointing  in  the 
jth  coordinate  direction.  Hence,  the  diagonal  entries  of  the  stress  tensor  a  refer  to  normal 
stresses,  e.  g.  due  to  compression,  and  the  off-diagonal  entries  represent  tangential  (shear) 
stresses.  The  Cauchy  stress  law  states  that  due  to  the  preservation  of  angular  momentum 
the  stress  tensor  a  is  symmetric  [13]. 

In  a  Newtonian  fluid  the  stress  tensor  is  assumed  to  depend  linearly  on  the  gradient  Vv 
of  the  velocity  v.  In  case  of  a  rigid  body  motion  the  stress  vanishes.  A  rotational  component 
of  the  local  motion  is  generated  by  the  antisymmetric  part  ^{T>v  —  (Vv)'^)  of  the  velocity 
gradient  Vv  :=  and  it  has  the  local  rotation  axis  V  x  u  and  local  angular 

velocity  |V  x  u|  [40].  Hence,  as  rotations  are  rigid  body  motions,  the  stress  only  depends  on 
the  symmetric  part  e[v]  :=  \{Vv-\-{T>vY')  of  the  velocity  gradient.  If  we  separate  compressive 
stresses,  reflected  by  the  trace  of  the  velocity  gradient,  from  shear  stresses  depending  solely 
on  the  trace- free  part  of  the  velocity  gradient,  we  obtain  the  constitutive  relation  of  an 
isotropic  Newtonian  fluid. 


^ij  h  (^shear)jj  T  A^c  (^bulk)jj  •  h 


dxi 


(2) 


where  /r  is  the  viscosity,  Kc  is  the  modulus  of  compression,  and  5ij  is  the  Kronecker  symbol. 

The  following  simple  configuration  serves  for  illustration.  We  consider  a  fluid  volume 
in  enclosed  between  two  parallel  plates  at  height  0  and  H,  where  the  vertical  direction 
normal  to  the  two  plates  points  along  the  a^^-coordinate  (cf.  Fig.  4).  Let  us  assume  the  lower 
plate  to  be  fixed  and  the  upper  plate  to  move  horizontally  at  speed  =  (uf ,  •  •  •  ,  0). 

Then,  the  velocity  field  v{x)  =  is  a  motion  field  consistent  with  the  boundary  conditions, 
and  the  resulting  stress  is  the  pure  shear  stress  acting  on  all  area  elements  parallel  to 
the  two  planes. 

Introducing  X  :=  Kc  —  and  denoting  the  jth  entry  of  the  Ah  row  of  e  by  Cjj,  one  can 
rewrite  (2)  as 

^ij  X5ij  ^  ]  c/jj;  -|-  , 

k 

or  in  matrix  notation  a  =  Atr(e)l  -|-  2/ie,  where  1  is  the  identity  matrix  and  e  =  e[v].  The 
parameter  A  is  denoted  Lame’s  first  coefficient.  The  local  rate  of  viscous  dissipation — the 
rate  at  which  mechanical  energy  is  locally  converted  into  heat  due  to  friction — can  now  be 
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computed  as 


diss[u] 


A 

2 

A 

2 


(tre[u])^  +  iiiv{e[vf) 

V  ^ 

\  i=l  /  i,j=l 


(Vjj  + 

4 


(3) 


where  we  abbreviated  Vij  =  ^-  To  see  this,  note  that  by  its  mechanical  dehnition,  the 
stress  tensor  a  is  the  hrst  variation  of  the  local  dissipation  rate  with  respect  to  the  velocity 
gradient,  i.  e.  a  =  ^/j^diss .  Indeed,  by  a  straightforward  computation  we  obtain 


Si^Dv)ij =  Atre(5p-  +  2/rep-  =  crp  . 

If  each  point  of  the  object  0{t)  at  time  t  G  [0, 1]  moves  at  the  velocity  v{x,t)  so  that  the 
total  deformation  of  0(0)  into  0{t)  can  be  obtained  by  integrating  the  velocity  field  v  in 
time,  then  the  accumulated  global  dissipation  of  the  motion  field  v  in  the  time  interval  [0, 1] 
takes  the  form 


Diss 


(n(t),C>(t))j^[op] 


0{t) 


diss[n]  da:  dt . 


(4) 


Here  tr(e[n]^)  measures  the  averaged  local  change  of  length  and  (tre[n])^  the  local  change  of 
volume  induced  by  the  transport.  Obviously  divn  =  tr(e[n])  =  0  characterizes  an  incom¬ 
pressible  fluid. 

Unlike  in  elasticity  models  (where  the  forces  on  the  material  depend  on  the  original 
configuration)  or  plasticity  models  (where  the  forces  depend  on  the  history  of  the  flow), 
in  the  Newtonian  model  of  viscous  fluids  the  rate  of  dissipation  and  the  induced  stresses 
solely  depend  on  the  gradient  of  the  motion  field  v  in  the  above  fashion.  Even  though  the 
dissipation  functional  (4)  looks  like  the  deformation  energy  from  linearized  elasticity,  if  the 
velocity  is  replaced  by  the  displacement,  the  underlying  physics  is  only  related  in  the  sense 
that  an  infinitisimal  displacement  in  the  fluid  leads  to  stresses  caused  by  viscous  friction, 
and  these  stresses  are  immediately  absorbed  via  dissipation,  which  reflects  a  local  heating. 

In  this  paper  we  address  the  problem  of  computing  geodesic  paths  and  distances  between 
non-rigid  shapes.  Shapes  will  be  modeled  as  the  boundary  contour  of  a  physical  object  that 
is  made  of  a  viscous  fluid.  The  fluid  flows  according  to  a  motion  field  v,  where  there  is  no  flow 
outside  the  object  boundary.  The  external  forces  which  induce  the  flow  can  be  thought  of 
as  originating  from  the  dissimilarity  between  consecutive  shapes.  The  resulting  Riemannian 
metric  on  the  shape  space,  which  defines  the  distance  between  shapes,  will  then  be  identified 
with  the  rate  of  dissipation,  representing  the  rate  at  which  mechanical  energy  is  converted 
into  heat  due  to  the  fluid  friction  whenever  a  shape  is  deformed  into  another  one. 


1.2  Related  work  on  shape  distances  and  geodesics 

Conceptually,  in  the  last  decade,  the  distance  between  shapes  has  been  extensively  studied 
on  the  basis  of  a  general  framework  of  the  space  of  shapes  and  its  intrinsic  structure.  The 
notion  of  a  shape  space  has  been  introduced  already  in  1984  by  Kendall  [25].  We  will  now 
discuss  related  work  on  measuring  distances  between  shapes  and  geodesics  in  shape  space. 
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particularly  emphasizing  the  relation  to  the  above  concepts  from  continuum  mechanics. 

An  isometrically  invariant  distance  measure  between  two  objects  Sa  and  Sb  in  (different) 
metric  spaces  is  the  Gromov-Hausdorff  distance  [23] ,  which  is  (in  a  simplified  form)  defined 
as  the  minimizer  of  |  \d{xi,  X2)  —  d(r/i,  r/2)|  over  all  maps  <p  ■.  Sa  ^  Sb  and 

■.  Sb  ^  Sai  matching  point  pairs  (xi,  X2)  in  <5^  with  pairs  (yi,  1/2)  in  Sb-  It  evaluates — 
globally  and  based  on  an  L°°-type  functional — the  lack  of  isometry  between  two  different 
shapes.  Memoli  and  Sapiro  [31]  introduced  this  concept  into  the  shape  analysis  community 
and  discussed  efficient  numerical  algorithms  based  on  a  robust  notion  of  intrinsic  distances 
d(-,  •)  on  shapes  given  by  point  clouds.  Bronstein  et  al.  incorporate  the  Gromov-Hausdorff 
distance  concept  in  various  classification  and  modeling  approaches  in  geometry  processing  [7] . 

In  [30]  Manay  et  al.  define  shape  distances  via  integral  invariants  of  shapes  and  demon¬ 
strate  the  robustness  of  this  approach  with  respect  to  noise. 

Gharpiat  et  al.  [10]  discuss  shape  averaging  and  shape  statistics  based  on  the  notion  of 
the  Hausdorff  distance  and  on  the  H^-noxm  of  the  difference  of  the  signed  distance  functions 
of  shapes.  They  study  gradient  flows  for  energies  defined  as  functions  over  these  distances 
for  the  warping  between  two  shapes.  As  the  underlying  metric  they  use  a  weighted 
metric,  which  weights  translational,  rotational,  and  scale  components  differently  from  the 
component  in  the  orthogonal  complement  of  all  these  transforms.  The  approach  by  Eckstein 
et  al.  [19]  is  conceptually  related.  They  consider  a  regularized  geometric  gradient  flow  for 
the  warping  of  surfaces. 

When  warping  objects  bounded  by  shapes  in  a  shape  tube  in  is  formed.  Delfour 
and  Zolesio  [15]  rigorously  develop  the  notion  of  a  Gourant  metric  in  this  context.  A  further 
generalization  to  classes  of  non-smooth  shapes  and  the  derivation  of  the  Euler-Lagrange 
equations  for  a  geodesic  in  terms  of  a  shortest  shape  tube  is  investigated  by  Zolesio  in  [48]. 

There  is  a  variety  of  approaches  which  consider  shape  space  as  an  infinite-dimensional 
Riemannian  manifold.  Michor  and  Mumford  [32]  gave  a  corresponding  definition  exempli¬ 
fied  in  the  case  of  planar  curves.  Yezzi  and  Mennucci  [43]  investigated  the  problem  that 
a  standard  L^-metric  on  the  space  of  curves  leads  to  a  trivial  geometric  structure.  They 
showed  how  this  problem  can  be  resolved  taking  into  account  the  conformal  factor  in  the 
metric.  In  [33]  Michor  et  al.  discuss  a  specific  metric  on  planar  curves,  for  which  geodesics 
can  be  described  explicitly.  In  particular,  they  demonstrate  that  the  sectional  curvature  on 
the  underlying  shape  space  is  bounded  from  below  by  zero  which  points  out  a  close  relation 
to  conjugate  points  in  shape  space  and  thus  to  only  locally  shortest  geodesics.  Younes  [44] 
considered  a  left-invariant  Riemannian  distance  between  planar  curves.  Miller  and  Younes 
generalized  this  concept  to  the  space  of  images  [34].  Klassen  and  Srivastava  [27]  proposed 
a  framework  for  geodesics  in  the  space  of  arclength  parametrized  curves  and  suggested  a 
shooting-type  algorithm  for  the  computation  whereas  Schmidt  et  al.  [37]  presented  an  alter¬ 
native  variational  approach. 

Dupuis  et  al.  [18]  and  Miller  et  al.  [35]  defined  the  distance  between  shapes  based  on  a 
flow  formulation  in  the  embedding  space.  They  exploited  the  fact  that  in  case  of  sufficient 
Sobelev  regularity  for  the  motion  field  v  on  the  whole  surrounding  domain  D,  the  induced 
flow  consists  of  a  family  of  diffeomorphisms.  This  regularity  is  ensured  by  a  functional 
Jg  Lv  ■  ndxdt,  where  L  is  a  higher  order  elliptic  operator  [39,44].  Thus,  if  one  considers 
the  computational  domain  Q,  to  contain  a  homogeneous  isotropic  fluid,  then  Lv  •  v  plays  the 
role  of  the  local  rate  of  dissipation  in  a  multipolar  fluid  model  [36] ,  which  is  characterized  by 
the  fact  that  the  stresses  depend  on  higher  spatial  derivatives  of  the  velocity.  Geometrically, 
JqLv  •  ndx  is  the  underlying  Riemannian  metric.  If  L  acts  only  on  e[v\  and  is  symmetric. 
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then  following  the  argnments  in  Section  1.1,  rigid  body  motion  invariance  is  incorporated 
in  this  multipolar  fluid  model.  Different  from  this  approach  we  conceptually  measure  the 
rate  of  dissipation  only  on  the  evolving  object  domain,  and  our  model  relies  on  classical 
(monopolar)  material  laws  from  fluid  mechanics  not  involving  higher  order  elliptic  operators. 
Under  sufficient  smoothness  assumptions  Beg  et  al.  derived  the  Euler-Lagrange  equations 
for  the  diffeomorphic  flow  field  in  [4].  To  compute  geodesics  between  hypersurfaces  in  the 
flow  of  diffeomorphism  framework,  a  penalty  functional  measures  the  distance  between  the 
transported  initial  shape  and  the  given  end  shape.  Vaillant  and  Glaunes  [41]  identified 
hypersurfaces  with  naturally  associated  two  forms  and  used  the  Hilbert  space  structures 
on  the  space  of  these  forms  to  define  a  mismatch  functional.  The  case  of  planar  curves  is 
investigated  under  the  same  perspective  by  Glaunes  et  al.  in  [22].  To  enable  the  statistical 
analysis  of  shape  structures,  parallel  transport  along  geodesics  is  proposed  by  Younes  et 
al.  [45]  as  the  suitable  tool  to  transfer  structural  information  from  subject-dependent  shape 
representations  to  a  single  template  shape. 

In  most  applications,  shapes  are  boundary  contours  of  physical  objects.  Fletcher  and 
Whitaker  [20]  adopt  this  view  point  to  develop  a  model  for  geodesics  in  shape  space  which 
avoids  overfolding.  Fuchs  et  al.  [21]  propose  a  Riemannian  metric  on  a  space  of  shape 
contours  motivated  by  linearized  elasticity,  leading  to  the  same  quadratic  form  (1)  as  in 
our  approach,  which  is  in  their  case  directly  evaluated  on  a  displacement  field  between  two 
consecutive  objects  from  a  discrete  object  path.  They  use  a  B-spline  parametrization  of 
the  shape  contour  together  with  a  finite  element  approximation  for  the  displacements  on 
a  triangulation  of  one  of  the  two  objects,  which  is  transported  along  the  path.  Due  to 
the  built-in  linearization  already  in  the  time-discrete  problem  this  approach  is  not  strictly 
rigid  body  motion  invariant,  and  interior  self-penetration  might  occur.  Furthermore,  the 
explicitly  parametrized  shapes  on  a  geodesic  path  share  the  same  topology,  and  contrary  to 
our  approach  a  cascadic  relaxation  method  is  not  considered. 

A  Riemannian  metric  in  the  space  of  3D  surface  triangulations  of  fixed  mesh  topology 
has  been  investigated  by  Kilian  et  al.  [26].  They  use  an  inner  product  on  time-discrete 
displacement  fields  to  measure  the  local  distance  from  a  rigid  body  motion.  These  local 
defect  measures  can  be  considered  as  a  geometrically  discrete  rate  of  dissipation.  Mainly 
tangential  displacements  are  taken  into  account  in  this  model.  Spatially  discrete  and  in  the 
limit  time-continuous  geodesic  paths  are  computed  in  the  space  of  discrete  surfaces  with  a 
fixed  underlying  simplicial  complex.  Recently,  Liu  et  al.  [28]  used  a  discrete  exterior  calculus 
approach  on  simplicial  complexes  to  compute  geodesics  and  geodesic  distances  in  the  space 
of  triangulated  shapes,  in  particular  taking  care  of  higher  genus  surfaces. 

1.3  Key  contributions 

The  main  contributions  of  our  approach  are  the  following: 

•  A  direct  connection  between  physics-motivated  and  geometry-motivated  shape  spaces 
is  provided,  and  an  intuitive  physical  interpretation  is  given  based  on  the  notion  of 
viscous  dissipation. 

•  The  approach  mathematically  links  a  pairwise  matching  of  consecutive  shapes  and 
a  viscous  flow  perspective  for  shapes  being  boundary  contours  of  objects  which  are 
represented  by  possibly  multi-labeled  images.  The  time  discretization  of  a  geodesic 


path  based  on  this  pairwise  matching  ensures  rigid  body  motion  invariance  and  a  1-1 
mapping  property. 

•  The  implicit  treatment  of  shapes  via  level  sets  allows  for  topological  transitions  and 
enables  the  computation  of  geodesics  in  the  context  of  partial  occlusion.  Robustness 
and  effectiveness  of  the  developed  algorithm  are  ensured  via  a  cascadic  multi-scale 
relaxation  strategy. 

2  The  variational  formulation 

Within  this  section,  in  2.1  we  put  forward  a  model  of  discrete  geodesics  as  a  finite  number 
of  shapes  Sk,  k  =  0, . . . ,  it',  connected  by  deformations  0^  :  Ok-i  which  are  optimal 

in  a  variational  sense  and  fulfill  the  hard  constraint  (pk{Sk-i)  =  Sk-  Subsequently,  in  2.3 
we  relax  this  constraint  using  a  penalty  formulation.  Afterwards,  based  on  a  viscous  fluid 
formulation,  in  2.4  we  introduce  a  model  for  geodesics  that  are  continuous  in  time,  and  in 
2.5  we  finally  show  that  the  latter  model  is  obtained  from  the  time-discrete  model  in  the 
limit  for  vanishing  time  step  size. 

2.1  The  time-discrete  geodesic  model 

As  already  outlined  above  we  do  not  consider  a  purely  geometric  notion  of  shapes  as  curves 
in  2D  or  surfaces  in  3D.  In  fact,  motivated  by  physics,  we  consider  shapes  S  as  boundaries 
do  of  sufficiently  regular,  open  object  domains  C>  C  for  d  =  2, 3.  Let  us  denote  by  S  a 
suitable  admissible  set  of  such  shapes  -  the  actual  shape  space.  Later,  in  Section  4.2,  this 
set  will  be  generalized  for  shapes  in  the  context  of  multi- labeled  images. 

Given  two  shapes  Sa,  Sb  in  S,  we  define  a  discrete  path  of  shapes  as  a  sequence  of  shapes 
(So,  (Si, ... ,  Sk  C  S  with  (So  =  Sa  and  Sk  =  Sb-  For  the  time  step  t  =  ^  the  shape  Sk 
is  supposed  to  be  an  approximation  of  S(tk)  for  tk  =  kr,  where  (‘S(t))jgjQ  is  a  continuous 
path  connecting  Sa  =  ‘S(O)  and  Sb  =  >5(1). 

Now,  we  consider  a  matching  deformation  :  Ok-i  for  each  pair  of  consecutive 

shapes  Sk-i  and  Sk  in  a  suitable  admissible  space  of  orientation  preserving  deformations 
D[C>fc_i]  and  impose  the  constraint  (pk{Sk-i)  =  Sk-  With  each  deformation  we  associate 
a  deformation  energy 


=  /  IF(D0fc)dX,  (5) 

JOk-1 

where  W  is  an  energy  density  which,  if  appropriately  chosen,  will  ensure  sufficient  regularity 
and  a  1-1  matching  property  for  a  deformation  0^  minimizing  ^Ideform  over  D[(9fc_i]  under  the 
above  constraint.  Analogously  to  the  axiom  of  elasticity,  the  energy  is  assumed  to  depend 
only  on  the  local  deformation,  reflected  by  the  Jacobian  T>(f)  :=  Yet,  different 

from  elasticity,  we  suppose  the  material  to  relax  instantaneously  so  that  object  Ok  is  again  in 
a  stress-free  configuration  when  applying  (/)k+i  at  the  next  time  step.  Let  us  also  emphasize 
that  the  stored  energy  does  not  depend  on  the  deformation  history  as  in  most  plasticity 
models  in  engineering. 

Given  a  discrete  path,  we  can  ask  for  a  suitable  measure  of  the  time-discrete  dissipation 
accumulated  along  the  path.  Here,  we  identify  this  dissipation  with  a  scaled  sum  of  the 
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accumulated  deformation  energies  ^^deforml^fc, along  the  path.  Furthermore,  the  inter¬ 
pretation  of  the  dissipation  rate  as  a  Riemannian  metric  motivates  a  corresponding  notion 
of  an  approximate  length  for  any  discrete  path.  This  leads  to  the  following  definition: 

Definition  1  (Discrete  dissipation  and  discrete  path  length).  Given  a  discrete  path  Sq, 
Si,  Sk  G  S,  the  total  dissipation  along  a  path  can  be  eomputed  as 


K 


D1SS7-(<Sq,  ,  .  .  .  ,  ^  ^  ^de/orm [^A:  1  — l]  •) 

L - /  j 


k=l 


where  4>k  is  a  minimizer  of  the  deformation  energy  ,  Sk-i]  over  Yy[Ok-\\  under  the 

constraint  fkiSk-i)  =  Sk-  Furthermore,  the  discrete  path  length  is  defined  as 


K 

Lt{So,  Si,  ,  Sk)  ■=  ^  V Sdeform[(l>k,Sk-l\  - 

k=l 

Let  us  make  a  brief  remark  on  the  proper  scaling  factor  for  the  time-discrete  dissipation. 
Indeed,  the  energy  S^,,form[<f>k,Sk-i]  is  expected  to  scale  like  r^.  Hence,  the  factor  ^  ensures 
a  dissipation  measure  which  is  conceptually  independent  of  the  time  step  size.  The  same 
holds  for  the  discrete  length  measure  ^/Sdetovm [0fc,  ‘5fc_i],  which  already  scales  like  r.  Thus 
Lr{So,  Si, ... ,  Sk)  indeed  reflects  a  path  length.  To  ensure  that  the  above-defined  dissipa¬ 
tion  and  length  of  discrete  paths  in  shape  space  are  well-defined,  a  minimizing  deformation 
fk  of  the  elastic  energy  Tdeform[‘, has  to  exist.  In  fact,  this  holds  for  objects  Ok-i  and 
Ok  with  Lipschitz  boundaries  Sk-i  and  Sk  for  which  there  exists  at  least  one  bi-Lipschitz 
deformation  fik  from  Ok-i  to  (9^  for  A:  =  1, . . . ,  iF  (i.  e.  fik  is  Lipschitz  and  injective  and  has 
a  Lipschitz  inverse).  The  associated  class  of  admissible  deformations  will  essentially  consist 
of  those  deformations  with  finite  energy.  Here,  we  postpone  this  discussion  until  the  energy 
density  of  the  deformation  energy  is  fully  introduced. 

With  the  notion  of  dissipation  at  hand  we  can  define  a  discrete  geodesic  path  following  the 
standard  paradigms  in  differential  geometry: 

Definition  2  (Discrete  geodesic  path).  A  discrete  path  Sq,  Si,  ,  Sk  in  a  set  of  admissible 
shapes  S  connecting  two  shapes  Sa  and  Sb  in  S  is  a  discrete  geodesic  if  there  exists  an 
associated  family  of  deformations  {(l)k)k=i,...,K  withfk  ^  F)[Ok-i\  and  (f)k{Sk-i)  =  Sk  such  that 
{(t>k,Sk)k=i,...,K  minimize  the  total  energy  'f2k=i^<ieform[fk,Sk-i]  over  all  intermediate  shapes 
Si, . . . ,  Sk-1  e  S  and  all  possible  matching  deformations  fi, . . . ,  ^k  with  fk  £  F)[Ok-i], 
Sk-i  =  dOk-i,  and  fik{Sk-i)  =  Sk  for  k  =  I, . . .  ,K . 

In  the  following,  we  will  inspect  an  appropriate  model  for  the  deformation  energy  density 
IF.  As  a  fundamental  requirement  for  the  time  discretization  we  postulate  the  invariance  of 
the  deformation  energy  with  respect  to  rigid  body  motions,  i.  e. 


^deform  [Q  °  +  ^)  ^k-l]  —  ^deform  [0fc  )  (6) 

for  any  orthogonal  matrix  Q  G  SO{d)  and  b  £  (the  axiom  of  frame  indifference  in  con¬ 
tinuum  mechanics).  From  this  one  deduces  that  the  energy  density  only  depends  on  the 
right  Cauchy-Green  deformation  tensor  'DcfF'Df,  i.  e.  there  is  a  function  IF  :  — >  M  such 

that  the  energy  density  IF  satisfies  IF(F)  =  W{F'^F)  for  all  F  £  Indeed,  if  (6)  holds 
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for  arbitrary  Sk-i,  4>k,  and  Q  E  SO{d),  then  we  have  to  have  W{QF)  =  W{F)  for  any 
Q  G  SO{d)  and  any  orientation  preserving  matrix  F  G  (in  particular,  F  =  T>(f)k{x)  for 
any  x  G  Ok-i).  By  the  polar  decomposition  theorem,  we  can  decompose  such  an  F  into 
the  product  of  an  orthogonal  matrix  Q  G  SO{d)  and  a  symmetric  positive  definite  matrix  C 

with  C  =  VF^F  and  Q  =  FVF^F~\  Thus,  W{F)  =  W{QVF^F)  =  W{VF^F)  so  that 
W{F)  can  indeed  be  rewritten  as  W{F^F),  where  W{C)  :=  W{'/C)  for  positive  definite 
matrices  C  G 

The  Cauchy-Green  deformation  tensor  geometrically  represents  the  metric  measuring  the 
deformed  length  in  the  undeformed  reference  configuration. 

For  an  isotropic  material  and  for  d  =  3  the  energy  density  can  be  further  rewritten  as  a  func¬ 
tion  W{Ii,  I2, 13)  solely  depending  on  the  principal  invariants  of  the  Cauchy-Green  tensor, 
namely  Ii  =  tr  (I) 0^1)0),  controlling  the  local  average  change  of  length,  I2  =  tr(cof(I)(/)^I)0)) 
(cofA  :=  det  AA~^),  reflecting  the  local  average  change  of  area,  and  I3  =  det(D0^D^), 
which  controls  the  local  change  of  volume.  For  a  detailed  discussion  we  refer  to  [14, 40] .  Let 
us  remark  that  tr(74"'"74)  coincides  with  the  Frobenius  norm  of  the  matrix  A  G  and 
the  corresponding  inner  product  on  matrices  is  given  hy  A  ■.  B  =  tx{AAB).  Furthermore,  let 
us  assume  that  the  energy  density  is  a  convex  function  of  !)(/),  coiVcj),  and  det  and  that 
isometries,  i.  e.  deformations  with  T>(f^ [x)T>(p{x)  =  1,  are  global  minimizers  [14].  For  the 
impact  of  this  assumption  on  the  time  discrete  geodesic  application  we  refer  in  particular  to 
the  second  row  in  Fig.  5,  which  provides  an  example  of  striking  global  isometry  preservation 
and  an  only  local  lack  of  isometry.  We  may  further  assume  kF(l)  =  W{d,d,  1)  =  0  without 
any  restriction.  An  example  of  this  class  of  energy  densities  is 

W{h,  I2,  h)  =  aj}  +  a2ll  +  T{h)  (7) 

with  p  >  1,  g  >  1,  tti  >  0,  q;2  >  0,  and  F  convex  with  r(/3)  ^  00  for  I3  ^  0  or  I3  ^  00, 
where  the  parameters  are  chosen  such  that  (hyh,  h)  =  {d,  d,  1)  is  the  global  minimizer  (cf. 
the  concrete  energy  density  defined  in  Appendix  A.l)  .  The  built-in  penalization  of  volume 

shrinkage,  i.  e.  W  00,  comes  along  with  a  local  injectivity  result  [3].  Thus,  the  sequence 
of  deformations  linking  objects  Ok-i  and  Ok  actually  represents  homeomorphisms  (which 
for  deformations  with  finite  energy  is  rigorously  proved  under  mild  assumptions  such  as 
sufficiently  large  p,  q,  certain  growth  conditions  on  F,  and  the  objects  embedded  in  a  very 
soft  instead  of  void  material  for  which  Dirichlet  boundary  conditions  are  prescribed).  We 
refer  to  [16],  where  a  similar  energy  has  been  used  in  the  context  of  morphological  image 
matching.  Let  us  remark  that  in  case  of  a  void  background,  self-contact  at  the  boundary 
is  still  possible  so  that  the  mapping  from  Sk-i  =  dOk-i  to  Sk  =  dOk  does  not  have  to  be 
homeomorphic.  With  the  interpretation  of  such  self-contact  as  a  closing  of  the  gap  between 
two  object  boundaries  in  the  sense  that  the  viscous  material  fiows  together,  our  model  allows 
for  topological  transitions  along  a  discrete  path  in  shape  space  [14]  (cf.  the  geodesic  from 
the  letter  A  to  the  letter  B  in  Fig.  1  for  an  example). 

2.2  An  existence  result  for  the  time-discrete  model 

Based  on  these  mechanical  preliminaries  we  can  now  state  an  existence  result  for  discrete 
geodesic  paths  for  a  suitable  choice  of  the  admissible  set  of  shapes  S  and  corresponding 
function  spaces  D[C>fc]  for  the  deformations  (pk,  k  =  I, . . . ,  K .  Note  that  the  known  regularity 
theory  in  nonlinear  elasticity  [3,12]  does  not  allow  to  control  the  Lipschitz  regularity  of  the 
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deformed  boundary  even  if  Sk-i  is  a  Lipschitz  boundary  of  the  elastic  domain  Ok-i- 

One  way  to  obtain  a  well-posed  formulation  of  the  whole  sequence  of  consecutive  variational 
problems  for  the  deformations  (pk  and  shapes  Sk  is  to  incorporate  the  required  regularity 
of  the  shapes  in  the  definition  of  the  shape  space.  Hence,  let  us  assume  that  S  consists  of 
shapes  S  which  are  boundary  contours  of  open,  bounded  sets  O  and  can  be  decomposed 
into  a  bounded  number  of  spline  surfaces  with  control  points  on  a  fixed  compact  domain. 
Furthermore,  the  shapes  are  supposed  to  fulfill  a  uniform  cone  condition,  i.  e.  each  point 
X  G  (S  is  the  tip  of  two  open  cones  with  fixed  opening  angle  a  >  0  and  height  r  >  0, 
one  contained  in  the  domain  O  and  the  other  in  the  complement  of  O.  On  such  object 
domains,  the  variational  problem  for  a  single  deformation  pk  connecting  shapes  Sk-i  and  Sk 
can  be  solved  based  on  the  direct  method  of  the  calculus  of  variations.  With  regard  to  the 
deformation  energy  integrand  in  (7),  the  natural  function  space  for  the  deformations  pk  is  a 
subset  of  the  Sobolev  space  [!]•  Let  us  take  into  account  an  explicit  function  F, 


namely  the  rational  function  F(/3)  =  «3  ^  j  ~  7-  Then,  in  d  =  3  dimensions,  for 

«!,  02,  as, /?,  7  >  0,  p,  g  >  3,  r  >  1  and  s  >  we  choose 


B[Ok-i]  :=  {P:Ok-i^R‘^\PeW^^P{Ok-i),coWpeL'^{Ok-i), 

detVp  E  L'^{Ok-i),detT>p  >  0  a.e.  in  Ok-i,  p{Ok-i)  =  Ok}  ■ 


Taking  into  account  this  space  of  admissible  deformations  for  each  k  G  {1, . . . ,  77}  leads  to 
a  well-defined  notion  of  dissipation  and  length  for  discrete  paths: 

Theorem  1  (Existence  of  a  discrete  geodesic).  Given  two  diffeomorphic  shapes  Sa  and  Sb 
in  the  above  shape  space  S,  there  exists  a  discrete  geodesic  So,  Si, ,  Sk  G  S  connecting 
Sa  and  Sb-  The  associated  deformations  pi, ... ,  pK  with  pk  G  'D[Ok-i\  for  k  =  1, . . .  ,K  are 
Holder  continuous  (that  is,  \p{x)  —  p{y)\  <\x  —  yp  for  some  7  G  (0, 1)  and  all  points  x,  y) 
and  locally  injective  in  the  sense  that  the  determinant  of  the  deformation  gradient  is  positive 
almost  everywhere. 

Proof:  To  prove  the  existence  of  a  discrete  geodesic  we  make  use  of  a  nowadays  classical 
result  from  the  vector- valued  calculus  of  variations.  Indeed,  applying  the  existence  results 
for  elastic  deformations  by  Ball  [2,3],  any  pair  of  consecutive  shapes  Sk-i  and  Sk  is  as¬ 
sociated  with  a  Holder  continuous  deformation  pk  G  D[C>fc_i]  with  detVpk  >  0  almost 
everywhere,  which  minimizes  the  deformation  energy  ^deform]',  *5^-1]  among  all  deformations 
p  G  D[Ofc_i].  Hence,  given  the  set  {pk)k=i,...,K  of  such  minimizing  deformations  for  fixed 
shapes  <Si, . . . ,  Sk,  we  can  compute  the  discrete  dissipation  }  ^deformi^fc,  along  the 

discrete  path  «Si, . . . ,  Sk- 

Now,  we  make  use  of  the  structural  assumption  on  the  shape  space  S.  The  space  of  all 
shapes  can  be  parametrized  with  finitely  many  parameters,  namely  the  control  points  of  the 
spline  segments.  These  control  points  lie  in  a  compact  set.  Also,  S  is  closed  with  respect  to 
the  convergence  of  this  set  of  parameters  since  the  cone  condition  is  preserved  in  the  limit 
for  a  convergent  sequence  of  spline  parameters. 

To  prove  that  a  minimizer  <Si, . . . ,  Sk  of  the  discrete  dissipation  Diss,-  exists,  we  first  observe 
that  DIsSt-  effectively  is  a  function  of  the  finite  set  of  spline  parameters.  Furthermore,  the 
set  of  admissible  spline  parameters  is  compact.  Hence,  it  is  sufficient  to  verify  that  Diss^ 
is  continuous.  For  this  purpose,  consider  shapes  Sk-i,  Sk  and  Sk-i,  Sk,  respectively.  Fur¬ 
thermore,  for  a  given  small  >  0  we  can  assume  the  spline  parameters  of  {Sk-i,Sk)  and 
{Sk-i,Sk)  to  be  close  enough  to  each  other  so  that  ioi  i  =  k  —  l,k  there  exists  a  bijective 
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deformation  ipi  :  Oi  ^  Oi  which  is  Lipschitz-continuous  and  has  a  Lipschitz-continuous 
inverse  ■0“^  with  —  l|ioo  +  “  ^lioo  —  for  a  d  <  ho-  Let  us  denote  by  (j),  0  the 

optimal  deformations  associated  with  the  dissipation  DissT-(«Sfc_i, <Sfc)  and  DissT-(«Sfc_i, «Sfc), 
respectively.  Using  the  optimality  of  0  and  defining  0  :=  o  0  o  ijjk-i  we  can  estimate 

JDisSr{Sk-i,Sk)-'DisSr{Sk-i,Sk)  =  ^  J  W{'D^)dx-^  j  W{'D(p)dx 

Ok-i  C>fc_i 

J  W{V(P)dx-^  J  W{V(P)dx 

Ok-i 

=  i  /  W{{Vi;-^o4>)V4>{V^Pk-ioi}k\))  \deiVi;-W  -  W{V(j))  da; . 

Ok-i 

Here,  we  have  applied  the  chain  rule  and  a  change  of  variables.  Taking  into  account  the 
explicit  form  of  the  integrand  and  the  above  assumption  on  ipk-i  and  V’fc,  we  can  estimate 
the  integrand  from  above  independently  of  S  by 

C'(ho)  +  |cofP^|^  +  |detD0|^  +  |(detD^)“^|*)  , 

where  C'(ho)  is  a  constant  solely  depending  on  ho.  Obviously,  this  pointwise  bound  itself  is 
integrable  for  (p  G  D((9fc-i).  Thus,  as  we  let  h  ^  0,  from  Lebesgue’s  theorem  we  deduce  that 

T>isSr{Sk-i,Sk)  -  'DisSr{Sk-i,Sk)  <  c(h) 

for  a  function  c  :  E"*"  ^  E  with  lim^^o  c(h)  =  0.  Exchanging  the  role  of  Sk-i,  Sk  and  Sk-i,  Sk 
we  obtain 


T>isSr{Sk-i,Sk)  -  T>isSr{Sk-i,Sk)  <  c{5) 

which  proves  the  required  continuity  of  the  dissipation  Diss,-.  Hence,  there  is  indeed  a  dis¬ 
crete  geodesic  <So, . . . ,  Sk-  D 


2.3  A  relaxed  formulation 

Computationally,  the  constraint  pkiSk-i)  =  Sk  for  a  1-1  matching  of  consecutive  shapes  is 
difficult  to  treat.  Furthermore,  the  constraint  is  not  robust  with  respect  to  noise.  Indeed, 
high  frequency  perturbations  of  the  input  shapes  Sa  and  Sb  might  require  high  deformation 
energies  in  order  to  map  Sa  onto  a  regular  intermediate  shape  or  to  obtain  <5^  as  the  image 
of  a  regular  intermediate  shape  in  a  1-1  manner.  Hence,  we  ask  for  a  relaxed  formulation 
which  allows  for  an  effective  numerical  implementation  and  is  robust  with  respect  to  noisy 
geometries.  At  first,  we  assume  that  the  complement  of  the  object  Ok-i  also  is  deformable, 
but  several  orders  of  magnitude  softer  than  the  object  itself.  Hence,  we  define 

St-.l  =  (<1  -  +  p  W(V4>t)  di  (8) 
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Figure  5:  Discrete  geodesic  for  two  different  examples  from  [21]  and  [11]  where  the  local 
rate  of  dissipation  is  color-coded  as  ■  ■.  In  the  bottom  example  the  local  preservation  of 

isometries  is  clearly  visible,  whereas  in  the  top  example  stretching  is  the  major  effect. 

for  deformations  (pk  now  defined  on  a  sufficiently  large  computational  domain  D.  For  simplic¬ 
ity  we  assume  0fc(x)  =  x  on  the  boundary  dVi.  This  renders  the  subproblem  of  computing 
an  optimal  elastic  deformation  well-posed  independent  of  the  current  shape.  For  5  =  0,  we 
obtain  the  original  model  and  suppose  that  at  least  a  sufficiently  smooth  extension  of  the 
deformation  on  a  neighborhood  of  the  shape  is  given. 

Now,  we  are  in  the  position  to  introduce  a  relaxed  formulation  of  the  pairwise  matching 
problem  by  adding  a  mismatch  penalty 

^match[0fc,‘5fc-l,Nfc]  =  VOl(C>fc_lA0^^(C)fc))  ,  (9) 

where  AAB  =  A\B(J B\A  defines  the  symmetric  difference  between  two  sets  and  vo^A)  = 
dx  is  the  d-dimensional  volume  of  the  set  A.  This  mismatch  penalty  replaces  the  hard 
matching  constraint  pki^k-i)  =  ^k-  Alternatively,  one  might  consider  the  mismatch  penalty 
vol((/)fc(C>fc_i)A(9fc),  but  as  we  will  see  in  Section  3.1,  the  form  (9)  is  computationally  more 
feasible  in  case  of  an  implicit  shape  description. 

Next,  in  practical  applications  shapes  are  frequently  defined  as  contours  in  images  and  usually 
not  given  in  explicit  parametrized  form.  Hence,  the  restriction  of  the  set  of  admissible  shapes 
to  piecewise  parametric  shapes,  which  we  have  taken  into  account  in  the  previous  section 
to  establish  an  existence  result  for  geodesic  paths,  is — from  a  computational  viewpoint — not 
very  appropriate  either.  If  we  allow  for  more  general  shapes  being  boundary  contours  of 
objects  in  images,  one  should  at  least  require  them  to  have  a  finite  perimeter.  Otherwise 
it  would  be  appropriate  to  decompose  the  initial  object  Oa  into  tiny  disconnected  pieces, 
shuffle  these  around  via  rigid  body  motions  (at  no  cost),  and  remerge  them  to  obtain  the 
final  object  Ob-  The  property  of  finite  perimeter  can  be  enforced  for  the  intermediate  shapes 
by  adding  the  object  perimeter  (generalized  surface  area  in  d  dimensions)  as  an  additional 
energy  term 

^area[‘5]  =  y  du  . 

Finally,  we  obtain  the  following  relaxed  definition  of  a  path  functional  for  a  family  of  defor¬ 
mations  and  shapes: 
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Figure  6:  Geodesic  paths  between  an  X  and  an  M,  without  a  contour  length  term  (i^  =  0,  top 
row) ,  allowing  for  crack  formation  (marked  by  the  arrows) ,  and  with  this  term  damping  down 
cracks  and  rounding  corners  (bottom  rows).  In  the  bottom  rows  we  additionally  enforced 
area  preservation  along  the  geodesic. 


Definition  3  (Relaxed  discrete  path  functional).  Given  a  sequence  of  shapes  {Sk)k=o,...,K  o,nd 
a  family  of  deformations  {(l>k)k=i,...,K  with  fk  '■  Ok-i  we  define  the  relaxed  dissipation 

as 

K 

£T[{4>k,Sk)k=l,...,K\  '■= 

k=l 

where  p,  u  are  parameters.  A  minimizer  of  this  energy  defines  a  relaxed  diserete  geodesic 
path  between  the  shapes  Sa  =  and  Sb  =  Sk- 

As  we  will  see  in  Section  2.5  below,  the  different  scaling  of  the  three  energy  components 
with  respect  to  the  time  step  size  r  will  ensure  a  meaningful  limit  for  r  ^  0. 

Fig.  6  shows  an  example  of  two  different  geodesics  between  the  letters  X  and  M,  demon¬ 
strating  the  impact  of  the  term  f^rea  controlling  the  {d  —  l)-dimensional  area  of  the  shapes. 


^deformW^kl  ij 


T 


+  V  6raatch[fik)  Gk-1  ,  Sk]  +  P  T  Sarea[Sk]  )  ,  (10) 


2.4  The  time- continuous  viscous  fluid  model 


In  this  section  we  discuss  geodesics  in  shape  space  from  a  Riemannian  perspective  and 
elaborate  on  the  relation  to  viscous  fluids.  This  prepares  the  identification  of  the  resulting 
model  as  the  limit  of  our  time  discrete  formulations  in  the  following  section.  A  Riemannian 
metric  ^  on  a  differential  manifold  Ad  is  a  bilinear  mapping  that  assigns  each  element  S  £  A4 
an  inner  product  on  variations  SS  oi  S.  The  associated  length  of  a  tangent  vector  SS  is  given 
by  ||5<S||  =  y^G{5S,  SS).  The  length  of  a  differentiable  curve  <S  :  [0, 1]  — Ad  is  then  defined 
by 


L[s]=  /  iiyi)iid(  = 


g(S(i),5(f))d(, 


where  S{t)  is  the  temporal  variation  of  S  at  time  t.  The  Riemannian  distance  between  two 
points  Sa  and  Sb  on  Ad  is  given  as  the  minimal  length  taken  over  all  curves  with  <S(0)  =  Sa 
and  <S(1)  =  Sb-  Hence,  the  shortest  such  curve  <S  :  [0, 1]  — >  Ad  is  the  minimizer  of  the  length 
functional  A  [5].  It  is  well-known  from  differential  geometry  that  it  is  at  the  same  time  a 
minimizer  of  the  cost  functional 


^(<S(t),<S(t))dt 
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and  describes  a  geodesic  between  Sa  and  Sb  of  minimum  length.  Let  us  emphasize  that 
a  general  geodesic  is  only  locally  the  shortest  curve.  In  particular  there  might  be  multiple 
geodesics  of  different  length  connecting  the  same  end  points. 

In  our  case  the  Riemannian  manifold  Ai  is  the  space  of  all  shapes  S  in  an  admissible  class  of 
shapes  S  (e.  g.  the  one  introduced  in  Section  2.1)  equipped  with  a  metric  Q  on  infinitesimal 
shape  variations.  As  already  pointed  out  above,  we  consider  shapes  S  as  boundary  contours 
of  deforming  objects  O.  Hence,  an  infinitesimal  normal  variation  5S  of  a  shape  S  =  dO  is 
associated  with  a  transport  field  v  :  O  ^  This  transport  field  is  obviously  not  unique. 
Indeed,  given  any  vector  field  w  on  O  with  w{x)  G  T^S  for  all  x  G  <S  =  dO  (where  T^S 
denotes  the  {d—  l)-dimensional  tangent  space  to  <S  at  a;),  the  transport  field  v  +  w  is  another 
possible  representation  of  the  shape  variation  5S.  Let  us  denote  by  V((5<S)  the  affine  space 
of  all  these  representations.  As  a  geometric  condition  for  v  G  V((5<S)  we  obtain  v  •  n[(S]  =  SS, 
where  n[<S]  denotes  the  outer  normal  of  S.  Given  all  possible  representations  we  are  interested 
in  the  optimal  transport,  i.  e.  the  transport  leading  to  the  least  dissipation.  Thus,  using  the 
definition  (1)  of  the  local  dissipation  rate  diss[n]  =  |  (tre[n])^  +  /rtr(e[n]^)  we  define  the 
metric  Q{SS,SS)  as  the  minimal  dissipation  on  motion  fields  v,  which  are  consistent  with 
the  variation  of  the  shape  SS: 

G{5S,SS)  :=  min  [  diss[n]dx=  min  /*  —  (trefnl)^  + /rtr(eM^)  da; .  (11) 

t;ev(<55)  Jq  vev(ss)  Jq  2 


Let  us  remark  that  we  distinguish  explicitly  between  the  metric  g{v,v)  :=  J^diss[n]da;  on 
motion  fields  and  the  metric  Q{SS,SS)  on  (the  different  space  of)  shape  variations,  which 
is  the  minimum  of  g{v,v)  over  all  motion  fields  consistent  with  SS.  Finally,  integration  in 
time  leads  to  the  total  dissipation 


min 

v(t)eV{S{t)) 


Diss 


(n(t),C>(t))jg[oi] 


g{s{t),s{t))dt 


to  be  invested  in  the  transport  along  a  path  (<S(t))te[o,i]  in  the  shape  space  S.  This  implies 
the  following  definition  of  a  time  continuous  geodesic  path  in  shape  shape: 

Definition  4  (Time-continuous  geodesic  path).  Given  two  shapes  Sa  and  Sb  in  a  shape 
spaee  S,  a  geodesic  path  between  Sa  and  Sb  is  a  curve  ((S(t))te[oq]  C  S  with  <S(0)  =  Sa  and 
«S(1)  =  Sb  which  is  a  local  solution  of 


min 

v(t)ev{S(t)) 


Diss 


{v{t),0{t)) 


te[o,i] 


among  all  differentiable  paths  in  S. 

Evidently,  one  has  to  minimize  over  all  motion  fields  v  in  space  and  time  which  are 
consistent  with  the  temporal  evolution  of  the  shape.  As  in  the  time-discrete  case,  we  can  relax 
this  property  and  consider  general  vector  fields  v  which  are  defined  at  time  t  on  the  domain 
0(t)  but  are  not  necessarily  consistent  with  the  evolving  shape.  The  lack  of  consistency  is 
instead  penalized  via  the  functional 


Top[(n(t),<S(t))t6[o,i]]  =  •  ?^[T‘5(t)]|  da. 


(12) 
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where  (1,  v{t))  is  the  underlying  space-time  motion  held  and  n[t,S{t)]  the  space-time  normal 
on  the  shape  tube  T  :=  Ute[oi](^’‘^(^))  ^  [0;  1]  H  denote  by  Xto  characteristic 

function  of  the  associated  (d-l-l)-dimensional  domain  tube  To  :=  Ute[o  i]  ^(^))  [0)  1] 

then — with  a  slight  misuse  of  notation — we  can  rewrite  this  functional  as 


£oF[{v{t),S{t))te[o,i]\  =  [  l^tXro  +  VxXto  •'^'I  dxdt.  (13) 

d(0,l)xR‘i 


Obviously,  there  is  a  similarity  to  TV-type  variational  approaches  in  optical  how  [5],  where 
V  is  the  optical  how  held  and  {t,x)  Xo(t)i^)  intensity  map  of  the  corresponding 

image  sequence. 

Additionally,  we  may  consider  a  further  regularization  term  on  the  tube  of  shapes,  which 
integrates  the  surface  area  E„^^[S{t)]  =  da  over  time  so  that  we  hnally  obtain  the 
time-continuous  path  functional 


T[(u(t),<S(t))te[o,i]]  =  /  [  diss[u]da;dt-hr7ToF[(t’(^),‘S(t))te[o,i]]- 

do  Jo(t) 


dadt.  (14) 


0  Js{t) 


Let  us  remark  that  the  second  and  the  third  energy  term  can  be  considered  as  anisotropic 
measures  of  area  on  the  space-time  tube  T.  Indeed,  the  last  term  integrates  the  (d  —  1)- 
dimensional  area  on  cross  sections  of  T  whereas  the  second  term  weights  the  area  element 
W {t,x)Xro  \  the  space  time  motion  held  (l,u). 


2.5  The  viscous  fluid  model  as  a  limit  for  r  ^  0 

We  now  investigate  the  relation  of  the  above-introduced  relaxed  discrete  geodesic  paths  and 
the  time  continuons  model  for  geodesics  in  shape  space.  For  this  purpose,  we  choose  the 
deformation  energy  in  such  a  way  that  the  Hessian  of  the  energy  Tdeform  with  respect  to  the 
deformation  of  an  object  O,  evaluated  at  the  identity  deformation  1,  coincides  up  to  a  factor 
I  with  the  dissipation  rate  or  metric  tensor  based  on  (1),  i.e. 

Hess£'deform[l,‘5](u,u)  =  2  /  diss[u]dx  (15) 

Jo 

for  any  velocity  held  v.  In  terms  of  the  energy  density  W  this  is  expressed  by  the  condition 

^^IF(1 -|- tA)|t=o  =  A(trA)^ -|- ^tr  ^(A -|- a"*")  ^  (16) 

for  the  second  derivative  of  W.  By  straightforward  computation  one  verihes  that  for  any 
local  dissipation  rate  (1)  one  can  hnd  a  nonlinear  energy  density  of  type  (7)  which  satishes 
(16).  This  is  detailed  in  Appendix  A.l,  expressing  the  free  parameters  of  the  deformation 
energy  density  (7)  in  terms  of  the  dissipation  parameters  A  and  p. 

Next,  let  us  introduce  the  following  notation.  Given  a  sequence  (So, . . .  ,Sk  of  shapes  and 
deformations  (pi, ...  ,4>k  with  pk  being  dehned  on  Ok-i,  we  introduce  a  temporally  piecewise 
constant  motion  held  and  a  time-continuous  deformation  held  p^  (which  interpolates 
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between  points  x  G  Ok-i  and  (t>k{x)  G  Ok)  by 


T 

4f;{t)  :=  (1  +  (t  - 

for  t  G  [tk-i,tk)  with  tk  =  kr.  The  corresponding  Eulerian  motion  field,  which  actually 
generates  the  flow,  is  then  given  by 


vAt)  :=  t.;  o  , 


Here,  we  assume  that  (j)^  is  injective. The  concatenation  with  its  inverse  is  only  needed  to 
obtain  the  proper  Eulerian  description  of  the  motion  field. 

For  decreasing  time  step  size  r,  we  are  interested  in  the  behavior  of  the  total  energy 
on  families  of  deformations  and  shapes,  given  by  the  time-discrete,  relaxed  model  from 
Definition  3,  and  its  relation  to  the  energy  £  on  motion  fields  and  shapes  in  space-time 
introduced  in  Definition  4.  In  fact,  if  we  evaluate  the  energy  on  a  family  of  deformations 
and  shapes,  where  the  deformations  are  induced  by  some  smooth  motion  field  v  and  the 
shapes  are  obtained  from  a  smooth  shape  tube  T  =  lJte[o  i](^’  ‘^(^))  regular  sampling,  we 
observe  convergence  to  the  time-continuous  energy  £  evaluated  on  v  and  T  as  postulated  in 
the  following  theorem: 


Theorem  2  (Limit  functional  for  vanishing  time  step  size).  Let  us  assume  that  (<S(t))te[o,i] 
is  a  smooth  family  of  shapes  and  consider  a  time  step  size  t  =  ^  with  K  ^  oo.  For  eaeh 
fixed  value  of  K  choose  Sk  =  S{kT)  for  k  =  0, ...  ,K.  Furthermore,  let  fi, ... ,  fx  be  a 
sequence  of  injective  deformations  with  fk  being  defined  on  Ok-i-  Finally,  assume  that  the 
associated  motion  field  Vr  converges  for  K  ^  oo  to  a  smooth  motion  field  v  on  the  space-time 
tube  UtG[o  i](^’ ^(^))-  Then  the  relaxed  diserete  path  functional  £r[{4>k,Sk)k=i,...,K]  converges 
to  the  time- continuous  path  functional  £’[(n(t),  »S(t))te[o,i]]  for  K  oo. 

We  conclude  that  our  variational  time  discretization  is  indeed  consistent  with  the  time- 
continuous  viscous  dissipation  model  of  geodesic  paths.  In  particular,  the  length  control 
based  on  the  first  invariant  Ii  of  Ffk  turns  into  the  control  of  infinitesimal  length  changes 
via  tr(e[n]^),  and  the  control  of  volume  changes  based  on  the  third  invariant  J3  of  Vfk  turns 
into  the  control  of  compression  via  tr(e[n])^  (cf.  Fig.  7  for  the  impact  of  these  two  terms  on 
the  shapes  along  a  geodesic  path).  Note  that  our  primal  interest  lies  in  the  case  1  since 
the  L^-type  optical  flow  term  is  supposed  to  just  act  as  a  penalty. 

Proof  (Theorem  2):  At  first,  let  us  investigate  the  convergence  behavior  of  the  sum  of 
deformation  energies  y  W[(9fc_i,  0^].  We  consider  a  second  order  Taylor  expansion 

around  the  identity  and  obtain 


W{Vfk)  =  IT(l)+TW^(l)(Dn^)  +  ylT,^A(l)(Dn^,Dn^)  +  0(r3) 

2  j2 

=  0  +  0  +  y^ir(l+i©J.J)|,.0  +  O(T!‘) 

=  +  ^tr  {(Vvt  +  (P«SP"))  +  0(e) 

=  r^diss[n:^]  -I-  O(r^) . 
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Figure  7:  Two  geodesic  paths  between  dumb  bell  shapes  varying  in  the  size  of  the  ends.  In 
the  top  example  the  ratio  \/ ji  between  the  parameters  of  the  dissipation  is  0.01  (leading 
to  rather  independent  compression  and  expansion  of  the  ends  since  the  associated  change 
of  volume  implies  relatively  low  dissipation),  and  100  in  the  bottom  example  (now  mass 
is  actually  transported  from  one  end  to  the  other).  The  underlying  texture  on  the  shape 
domains  Oq,  ■  ■  ■ ,  Ok-i  is  aligned  to  the  transport  direction,  and  the  absolute  value  of  the 
velocity  v  is  color-coded  as  ft  Ad. 


Here,  we  have  used  that  the  identity  deformation  is  the  minimizer  of  VF(-)  with  VF(1)  =  0 
as  well  as  the  relation  between  W  and  diss  from  (16).  Now,  summing  over  all  deformation 
energy  contributions  yields 


K 


lim  - 

K^oo  ^  T 
k=l 


W{Dct>k)dx 


dx  dt 


so  that  we  recover  the  viscous  dissipation  in  the  limit. 

Next,  we  investigate  the  limit  behavior  of  the  sum  of  mismatch  penalty  functionals  for 
vanishing  time  step  size.  In  a  neighborhood  of  the  shape  Sk-i,  let  us  for  x  G  S^-i  define  the 
local  and  signed  thickness  function  (cf.  Fig.  8) 


8k{x)  :=sup{s  :  (f)k{x sn[Sk-i]{x))  G  Ok) 

of  the  mismatch  set  A0^^(C>fc)  (recall  that  (pk  is  extended  outside  Ok-i).  Then,  we 
obtain 


vol{Ok-iA(p,^^{Ok))=  f  |4(a::)|do-Fo(T) .  (17) 

^  <^k-l 

Furthermore,  we  connect  the  shapes  Sk-i  =  S(tk-i)  and  Sk  =  S(tk)  via  a  ruled  surface 
T^i^ied;  For  x  G  Sk-i  we  suppose  a  vector  rk{x)  G  with  =  0{t)  to  be  defined  by  the 
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properties  rk{x)  _L  T^Sk-i  and  x  +  rk{x)  G  Sk-  Then  define 


J-ruled  ^ ;  f  g  [tk_l,tk],X  G  Sk-1^  ■ 

Obviously,  7^’^’^'®'^  approximates  the  continuous  tube  %  :=  up  to  terms  of 

the  order  0{t^).  We  denote  by  nk[tk-i,Sk-i\{x)  the  normal  vector  on  the  ruled  surface 
^ruied  point  X  G  Sk-i-  In  particular,  nk[tk-i,Sk-i]{x)  T  (0,t(;)  Vw  G  T^Sk-i  and 
nk[tk-i,  Sk-i]{x)  ±  {T,rk{x)).  From  these  properties  we  get  that 

|(T,rfc(a;)  -  5kix)n[Sk-i]{x))  ■  nk[tk-i,Sk-i]ix)\  =  r  |(l,n^(x))  •  nk[tk-i,Sk-i]{x)\  +  o(t)  .(18) 
Next,  by  an  elementary  geometric  argument  for 
Ikix)  :=  \/r2  +  |rfc(a;)P  , 

ek(x)  :=  {T,rk{x)  -  4(x)n[<Sfc_i](a:))  •  nk[tk-i,Sk-i]{x) 
we  obtain  that  and  hence 

T 

Using  this  relation  together  with  (18)  and  taking  into  account  further  standard  approxima¬ 
tion  arguments  we  obtain 


1(1,  n(x))  ■n[t,S{t)]{x)\da=  /  |(l,n^(x))  •  nfc[4-i, ‘5fc-i](a:)|  da -h  o(t) 


'Tfc 


/  '^ruled 


\{l,v';{x))  •  nk[tk-i.,Sk-i]{x)\lk{x)  da  o(t) 


'Su- 


'Sk- 


-\{T,rk{x)  -  5k{x)n[Sk-i\{x))  ■  nk[tk-i,Sk-i\{x)\lk{x)  da -h  o(t) 

T 

|4(a:)|  da  -I-  o(t) 


so  that  by  (17)  we  finally  arrive  at  the  desired  result 


vol(C)fc-iA^^UC)fc))  =  /  \{l,v{x))  ■n[t,S{t)]{x)\da  +  o{T) . 


Finally,  the  sum  of  shape  perimeters,  Xfc=i  '^^area[‘5fc],  obviously  converges  to  the  time  integral 
of  the  perimeters, 


□ 


da  dt 

lo  JS(t) 

so  that  we  have  verified  the  postulated  convergence. 


Let  us  remark  that  we  do  not  prove  F-convergence  of  the  relaxed  discrete  path  functional  as 
the  time  step  size  approaches  zero.  Here,  the  issue  of  compactness  of  the  family  of  shapes 
and  deformations  with  finite  energy  as  well  as  the  lower  semi-continuity  are  open  problems. 
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Figure  8:  Sketch  of  the  mismatch  between  shapes  and  motion  fields.  The  left  sketch  illus¬ 
trates  the  quantities  from  the  proof  for  a  geodesic  path  of  2D  shapes,  and  the  middle  shape 
shows  a  close-up.  The  right  graph  shows  the  corresponding  variables  in  space-time. 


Particularly  the  influence  of  the  anisotropic  area  measures  on  the  shape  tubes  in  space-time 
on  the  compactness  of  a  sequence  of  discrete  geodesics  for  vanishing  time  step  size  r  is  one 
of  the  major  challenges. 


3  The  numerical  algorithm 

In  this  section  we  deal  with  the  derivation  of  a  numerical  scheme  to  effectively  compute  the 
discrete  geodesic  paths.  In  Section  3.1  we  will  introduce  a  regularized  level  set  description  of 
shape  contours  and  rewrite  the  different  energy  contributions  of  (10)  in  terms  of  level  sets. 
Then,  a  spatial  finite  element  discretization  for  the  level  set-based  shape  description  and  the 
deformations  (pk  is  investigated  in  3.2.  Finally,  a  sketch  of  the  resulting  numerical  algorithm 
is  given  in  3.3. 

3.1  Regularized  level  set  approximation 

To  numerically  solve  the  minimization  problem  for  the  energy  (10),  we  assume  the  object 
domains  Ok  to  be  represented  by  zero  super  level  sets  {x  G  D  :  Uk{x)  >  0}  of  a  scalar  function 
Mfc  :  D  ^  E  on  a  computational  domain  D  C  E^.  Similar  representations  of  shapes  have  been 
used  for  shape  matching  and  warping  in  [10,24].  We  follow  the  approximation  proposed  by 
Chan  and  Vese  [9]  and  encode  the  partition  of  the  domain  D  into  object  and  background 
in  the  different  energy  terms  via  a  regularized  Heaviside  function  Hg{uk).  As  in  [9]  we 
consider  the  function  H^^x)  ■=  ^  +  ^  arctan  (|)  ,  where  e  is  a  scale  parameter  representing 
the  width  of  the  smeared-out  shape  contour.  Hence,  the  mismatch  energy  is  replaced  by  the 
approximation 


^Ltch[(Pk,  Uk-i,Uk]  =  /  {He{Uk  o  <pk)  -  He{Uk-i))  dx  , 

Jn 

and  the  area  of  the  fcth  shape  Sk  is  replaced  by  the  total  variation  of  o  Uk, 


C.N=  /  |vff,K)|dx 

Jn 


(19) 


(20) 
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In  the  expression  for  the  relaxed  elastic  energy  (8)  we  again  replace  the  characteristic  function 
Xok-i  Hs{uk)  and  obtain 

[  {il-5)H,{u,.,)  +  S)W{V^k)dx,  (21) 

Jci 

where  h  =  10“^  in  our  implementation.  Let  us  emphasize  that  in  the  energy  minimization 
algorithm,  the  guidance  of  the  initial  zero  level  lines  towards  the  final  shapes  relies  on  the 
nonlocal  support  of  the  derivative  of  the  regularized  Heaviside  function  (cf.  [8]).  Finally,  we 
end  up  with  the  approximation  of  the  total  energy, 

Sr^[{4>k,  Uk)k=l,...,K]  =  ^  Uk-l]  +  V^Ltch[4>k,  Uk-l,Uk]  +  ■  (22) 

k=l  ^ 

In  our  applications  we  have  chosen  values  for  rj  between  20  and  200  and  u  either  zero  or 
0.001  (except  for  Fig.  6,  where  i'  =  0.05).  Within  these  ranges,  the  shapes  along  the  discrete 
geodesics  are  relatively  independent  of  the  actual  parameter  values.  The  Lame  coefficients 
are  A  =  /r  =  1  apart  from  Fig.  7.  The  essential  formulas  for  the  variation  of  the  different 
energies  can  be  found  in  Appendix  A. 2. 

Note  that  in  order  to  be  a  proper  approximation  of  the  model  with  sharp  contours,  £ 
should  be  smaller  than  the  shape  variations  between  consecutive  shapes  along  the  discrete 
geodesic.  Only  in  that  case,  the  integrand  of  (19)  is  one  on  most  of  (9fc-iA0^^(Clfc).  Conse¬ 
quently,  as  T  ^  0,  £  has  to  approach  zero  at  least  at  the  same  rate. 

3.2  Finite  element  discretization  in  space 

For  the  spatial  discretization  of  the  energy  6^^  (^22)  the  finite  element  method  has  been 

applied.  The  level  set  functions  and  the  different  components  of  the  deformations  are 
represented  by  continuous,  piecewise  multilinear  (trilinear  in  3D  and  bilinear  in  2D)  finite 
element  functions  Uk  and  on  a  regular  grid  superimposed  on  the  domain  D  =  [0,  l]*^.  For 
the  ease  of  implementation  we  consider  dyadic  grid  resolutions  with  2^  -|-  1  vertices  in  each 
direction  and  a  grid  size  h  =  2“^.  In  2D  we  have  chosen  L  =  7, . . . ,  10  and  in  3D  L  =  7. 
Single  level  minimization  algorithm.  For  fixed  time  step  r  and  fixed  spatial  grid  size  h,  let  us 
denote  by  the  discrete  total  energy  depending  on  the  set  of  K  discrete 

deformations  d>i, . . . ,  d>x  and  A'  +  1  discrete  level  set  functions  Uq,  . . . ,  Uk,  where  Uq  and 
Uk  describe  the  shapes  Sa  and  Sb  and  are  fixed.  This  is  a  nonlinear  functional  both  in 
the  discrete  deformations  (due  to  the  concatenation  7/^  o  with  the  discrete  level  set 
function  Uk  and  the  nonlinear  integrand  W{-)  of  the  deformation  energy  f^deLm)  as  well  as  in 
the  discrete  level  set  functions  Uk  (due  to  the  concatenation  with  the  regularized  Heaviside 
function  i7e(-))-  In  onr  energy  relaxation  algorithm  for  fixed  time  step  and  grid  size,  we 
employ  a  gradient  descent  approach.  We  constantly  alternate  between  performing  a  single 
gradient  descent  step  for  all  deformations  and  one  for  all  level  set  functions.  The  step  sizes  are 
chosen  according  to  Armijo’s  rule.  If  the  actually  observed  energy  decay  in  one  step  is  smaller 
than  I  of  the  decay  estimated  from  the  derivative  (the  Armijo  condition  is  then  declared 
to  be  violated),  then  the  step  size  is  halved  for  the  next  trial,  else  it  is  doubled  as  often  as 
possible  without  violating  the  Armijo  condition.  This  simultaneous  relaxation  with  respect 
to  the  whole  set  of  discrete  deformations  and  discrete  level  set  functions  (representing  the 
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shapes),  respectively,  already  outperforms  a  simple  nonlinear  Gauss-Seidel  type  relaxation 
(cf.  Fig.  3).  Nevertheless,  the  capability  to  identify  a  shortest  path  between  complicated 
shapes  depends  on  an  effective  multi-scale  relaxation  strategy  (see  below). 

Numerical  quadrature.  Integral  evaluations  in  the  energy  descent  algorithm  are  performed  by 
Gaussian  quadrature  of  third  order  on  each  grid  cell.  For  various  terms  we  have  to  evaluate 
pullbacks  17  o  $  of  a  discretized  level  set  function  17  or  a  test  function  under  a  discretized 
deformation  <!>.  Let  us  emphasize  that  quadrature  based  on  nodal  interpolation  of  17  o  <I> 
would  lead  to  artificial  displacements  near  the  shape  edges  accompanied  by  strong  artificial 
tension.  Hence,  in  our  algorithm,  if  d>(x)  lies  inside  for  a  quadrature  point  x,  then  the 
pullback  is  evaluated  exactly  at  x.  Otherwise,  we  project  d>(a:)  back  onto  the  boundary  of 
O  and  evaluate  U  at  that  projection  point.  This  procedure  is  important  for  two  reasons: 
First,  if  we  only  integrated  in  regions  for  which  $(a;)  €  O,  we  would  induce  a  tendency  for 
$  to  shift  the  domain  outwards  until  $(G)  fl  O  =  0,  since  this  would  yield  zero  mismatch 
penalty.  Second,  for  a  gradient  descent  to  work  properly,  we  need  a  smooth  transition  of  the 
energy  if  a  quadrature  point  is  displaced  outside  fl  or  comes  back  in.  By  the  form  of  the 
mismatch  penalty,  this  implies  that  the  discrete  level  set  functions  Uk  have  to  be  extended 
continuously  outside  Id.  Backprojecting  'I>(x)  onto  the  boundary  just  emulates  a  constant 
extension  of  17^  perpendicular  to  the  boundary. 

Cascadic  multi-scale  algorithm.  The  variational  problem  considered  here  is  highly  nonlinear, 
and  for  fixed  time  step  size  the  proposed  scheme  is  expected  to  have  very  slow  convergence; 
also  it  might  end  up  in  some  nearby  local  minimum.  Here,  a  multi-level  approach  (initial 
optimization  on  a  coarse  scale  and  successive  refinement)  turns  out  to  be  indispensable  in 
order  to  accelerate  convergence  and  not  to  be  trapped  in  undesirable  local  minima.  Due  to 
our  assumption  of  a  dyadic  resolution  2^  -|-  1  in  each  grid  direction,  we  are  able  to  build  a 
hierarchy  of  grids  with  2^  -|- 1  nodes  in  each  direction  for  I  =  L, . . .  ,0.  Via  a  simple  restriction 
operation  we  project  every  finite  element  function  to  any  of  these  coarse  grid  spaces.  Starting 
the  optimization  on  a  coarse  grid,  the  results  from  coarse  scales  are  successively  prolongated 
onto  the  next  grid  level  for  a  refinement  of  the  solution  [6].  Hence,  the  construction  of  a  grid 
hierarchy  allows  to  solve  coarse  scale  problems  in  our  multi-scale  approach  on  coarse  grids. 
Since  the  width  £  of  the  diffusive  shape  representation  o  Uk  should  naturally  scale  with  the 
grid  width  h,  we  choose  e  =  h.  Likewise,  we  first  start  with  a  coarse  time  discretization  and 
successively  add  intermediate  shapes.  At  the  beginning  of  the  algorithm,  the  intermediate 
shapes  are  initialized  as  one  of  the  end  shapes. 

On  a  3  GHz  Pentium  4,  still  without  runtime  optimization,  2D  computations  for  L  =  8  and 
K  =  8  require  ~  1  h.  Based  on  a  parallelized  implementation  we  observed  almost  linear 
scaling. 

3.3  A  sketch  of  the  algorithm 

The  entire  algorithm  in  pseudo  code  notation  reads  as  follows  (where  bold  capitals  represent 
vectors  of  nodal  values  and  the  2^  +  1  shapes  on  time  level  j  are  labeled  with  the  superscript 

j)- 


EnergyRelaxation  (Lstart,  Ge„d)  { 
for  time  level  j  =  jo  to  J  { 

K  =  2^';  =  [/start;  Uj,  =  [/e„d 

if  (j  =  jo)  { 
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Figure  9:  Geodesic  path  between  a  cat  and  a  lion,  with  the  local  rate  of  dissipation  inside 
the  shapes  (So, ,  Sk-i  color-coded  as  ■  ■  (middle)  and  a  transparent  slicing  plane  with 

texture-coded  level  lines  of  the  level  set  representation  (bottom). 


initialize  =  1,  Uf  =  17^,  z  =  1, . . . ,  K 

}  else  { 

initialize  =  1  +  “  1),  ° 


UL  =  ur.  Ul,_,  =  ur  o  i  =  1, . . . , 


K. 

2  ’ 


} 

restrict  Uj,  for  alH  =  1, . . . ,  iF  onto  the  coarsest  grid  level  Iq; 
for  grid  level  I  =  Iq  to  L  { 
for  step  /c  =  0  to  fcmax  { 

perform  a  gradient  descent  step 


with  Armijo  step  size  control  for  r; 


{^i)i=i,...,K  =  -  r  grad(^oid) . $i)i=i,...,ic] 

1  for  r 

perform  a  gradient  descent  step 

=  (U°'^)i=i,„.,i^  -Tgrad(uoid),^^ 

with  Armijo  step  size  control  for  r; 

} 

if  (/  <  L)  prolongate  Uj,  for  alH  =  1, . . . ,  K  onto  the  next  grid  level; 


K 


4  Experimental  results  and  generalizations 

We  have  computed  discrete  geodesic  paths  for  2D  and  3D  shape  contours.  The  method  is 
both  robust  and  flexible  due  to  the  underlying  implicit  shape  description  via  level  sets,  cf. 
Fig.  1,  5,  7,  9,  10,  and  11.  Indeed,  neither  topologically  equivalent  meshes  on  the  end  shapes 
are  required,  nor  need  the  shapes  themselves  be  topologically  equivalent. 

In  what  follows  let  us  focus  on  a  number  of  different  applications  of  the  developed  com- 
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Figure  10:  Geodesic  path  between  the  hand  shapes  ni336  and  m324  from  the  Princeton 
Shape  Benchmark  [38].  Two  different  views  are  presented  in  the  first  two  rows.  The  bottom 
row  shows  the  local  dissipation  color-coded  on  slices  through  the  hand  shapes. 


putational  tool  and  suitable  extensions.  A  slight  modification  of  the  matching  condition, 
presented  in  Section  4.1,  will  allow  the  computation  of  discrete  geodesic  paths  in  case  of 
partial  occlusion  of  one  of  the  end  shapes.  Section  4.2  deals  with  the  fact  that  frequently, 
physical  objects  consists  of  different  regions.  Along  a  geodesic  path,  each  of  these  regions  has 
to  be  transported  consistently  from  one  object  onto  the  corresponding  region  in  the  other 
object.  Based  on  the  concept  of  multi-labeled  images  which  implicitly  represent  such  phys¬ 
ical  objects,  Section  4.2  generalizes  our  concept  of  geodesics  correspondingly.  Furthermore, 
the  computation  of  distances  between  groups  of  shapes  can  be  used  for  shape  statistics  and 
clustering,  which  will  be  considered  in  Section  4.3.  Finally,  we  will  show  in  Section  4.4  that 
already  for  simple  shapes  such  as  letters  there  might  be  multiple  (locally  shortest)  geodesics 
between  pairs  of  shapes.  The  shown  examples  will  not  only  give  some  deeper  insight  into 
the  structure  of  the  shape  space,  but  also  illustrate  the  stability  of  our  computational  results 
with  respect  to  geometric  shape  variations. 

4.1  Computing  geodesics  in  case  of  partial  occlusion 

In  many  shape  classification  applications,  one  would  like  to  evaluate  the  distance  of  a  partially 
occluded  shape  from  a  given  template  shape.  For  example  in  [17]  such  a  problem  has 
been  studied  in  the  context  of  joint  registration  of  multiple,  partially  occluded  shapes.  Our 
geodesic  model  can  be  adapted  to  allow  for  partial  occlusion  of  one  of  the  input  shapes.  Let 
us  suppose  that  the  domain  Oq  associated  with  the  shape  Sa  =  OOq  is  partically  occluded. 
Thus,  we  replace  the  first  term  in  the  sum  of  mismatch  penalty  functionals  by 

^match[01)  >50,  <Sl]  =  VOl(Clo  \  01  ^{Ol)) 
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Figure  11:  Top:  Real  video  sequence  of  a  white  blood  cell  (courtesy  Robert  A.  Freitas, 
Institute  for  Molecular  Manufacturing,  California,  USA).  Middle:  Discrete  geodesic  between 
the  corresponding  end  shapes.  Bottom:  Pushforward  of  the  first  image  under  a  concatenation 
of  the  deformations  connecting  consecutive  shapes  along  the  discrete  geodesic.  Note  that 
the  geodesic  interpolation  is  similar  to  the  actual  shape  deformation  observed  in  the  video. 

KKf: 


Figure  12:  A  discrete  geodesic  connecting  different  poses  of  a  matchstick  man  can  be  com¬ 
puted  (from  left  to  right  starting  with  the  second),  even  though  part  of  one  arm  and  one  leg 
of  So  (left)  are  occluded. 


and  do  not  penalize  areas  of  0]"^((Pi),  which  are  not  covered  by  the  (partially  occluded) 
domain  Oq-  Hence,  the  energy  will  favor  discrete  paths  in  shape  space  which  are  pairwise 
in  a  1-1  correspondence  except  for  the  very  first  pair,  where  only  an  approximate  inclusion 
of  Oq  in  is  intended.  For  this  purpose,  in  the  numerical  implementation  we  insert 

a  masking  function  i/£(Erosion£[Mo])  and  obtain 

“i]  =  /  o  0i)  -  He{uo)f  i7£(Erosion£-[Mo])  dx  . 

Ja 

Here  Erosion^  is  an  erosion  operator  acting  on  the  image  uq  and  eroding  the  domain  Oq  by 
a  width  i.  Eurthermore,  e  is  chosen  roughly  of  the  same  size  as  e  (we  actually  choose  e  =  e). 
This  modification  improves  the  robustness  of  the  descent  scheme  since  it  does  not  penalize 
deviations  of  the  pulled  back  level  set  function  Ui  o  (pi  from  Uq  in  the  interface  region  between 
the  occluded  and  non-occluded  parts  of  Oq.  An  application  of  this  modified  scheme  is  shown 
in  Fig.  12. 
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Figure  13:  Discrete  geodesic  between  the  straight  and  the  folded  bar  from  Fig.  2,  where 
the  black  region  of  the  initial  shape  in  the  top  row  is  constrained  to  be  matched  to  the 
black  region  of  the  final  shape.  The  bottom  row  shows  a  color-coding  of  the  corresponding 
viscous  dissipation.  Due  to  the  strong  difference  in  relative  position  of  the  black  region 
between  initial  and  end  shape,  the  intermediate  shapes  exhibit  a  strong  asymmetry  and  high 
dissipation  in  the  light  grey  region  near  both  ends  of  the  bar. 

4.2  Geodesics  between  multi- labeled  images 

Only  taking  into  account  shapes  which  are  outer  boundary  contours  S  =  dO  of  open  objects 
0  C  is  rather  limiting  in  some  applications.  While  the  contours  of  an  object  O a  are 
correctly  mapped  onto  the  contours  of  an  object  Ob  via  the  geodesic  between  Sa  =  OOa 
and  Sb  =  OOb,  the  viscous  fluid  model  imposes  no  restriction  on  the  path- generating  flow  in 
the  object  interior  (apart  from  the  property  that  it  should  minimize  the  viscous  dissipation). 
However,  one  might  often  want  certain  regions  of  one  object  Oa  to  be  mapped  onto  particular 
regions  in  another  object  Ob-  Generally  speaking,  real  world  shapes  or  objects  are  often 
characterized  as  a  composition  of  different  structures  or  components  with  a  particular  relative 
position  to  each  other.  A  geodesic  or  a  general  path  between  two  such  shapes  should  of  course 
match  corresponding  structures  with  each  other,  and  a  change  in  relative  position  of  these 
subcomponents  naturally  has  to  contribute  to  the  path  length. 

As  an  example,  let  us  reconsider  the  discrete  geodesic  between  the  straight  and  the 
folded  bar  in  Fig.  2.  The  initial  and  the  final  shape  contain  no  additional  information  about 
any  internal  structures  so  that  the  deformation  strength  and  the  induced  dissipation  along 
the  geodesic  path  are  distributed  evenly  over  the  whole  object,  in  particular  generating 
symmetric  intermediate  shapes.  However,  if  we  prescribe  the  original  and  the  final  location 
for  some  internal  region  of  the  bar,  the  dissipation-minimizing  flow  may  look  very  different  if 
the  additional  constraints  are  not  consistent  with  the  geodesic  flow  without  constraints  (cf. 
Fig.  13). 

For  these  reasons  we  would  like  to  extend  our  approach  to  allow  for  more  general  shapes 
that  may  be  composed  of  a  number  of  subcomponents.  Since  we  can  interpret  also  images 
as  collections  of  different  shapes  or  objects,  the  computation  of  geodesics  between  (multi- 
labeled)  images  nicely  fits  into  this  setting  as  well. 

The  extension  is  very  simple:  Instead  of  a  geodesic  between  just  two  shapes  Sa  =  OOa 
and  Sb  =  OOb,  we  now  seek  a  geodesic  path  {S^{t))i=i^,,,^rn  =  (5Cl*(^))«=i, ^  ^  [0)  !]>  be¬ 
tween  two  collections  of  shapes,  each  of  them  consisting  of  m  separate  shapes,  {SA)i=i,...,m  = 
{dOA)i=i,...,m  and  {S'^B)i=i,--,m  =  The  geodesic  path  is  supposed  to  be  gener¬ 
ated  by  a  joint  motion  field  v{t)  :  Ui=i  The  single  objects  0^{t)  can  then 

be  regarded  as  the  subcomponents  of  an  overall  object  Ui=i  m  The  total  dissipation 


27 


along  the  path  is  measured  exactly  as  before  by 


Diss[n]  =  f  f  ^  (tre[n])^  +  ptr(e[n]^)  dxdt . 

Jo  2 

This  naturally  translates  to  the  objective  functional  of  the  discrete  geodesic  with  K  +  1 
intermediate  shape  collections  k  =  0, . . . ,  K, 

K  ^  r 

k=l  k=l  Ui=l,...,Tn  ^k—1 

where  the  deformations  satisfy  the  constraints  <pk{Sl_i)  =  Si  for  k  =  i  = 

1, . . . , m,  and  S^  =  S\,  S}^  =  S^,  i  =  1, . . .  ,m. 

The  corresponding  relaxed  formulation  then  has  to  include  multiple  mismatch  penal¬ 
ties  (one  for  every  constraint),  and  as  before,  we  incorporate  a  regularization  of  the  shape 
perimeter  (or  generalized  surface  area  in  d  dimensions)  so  that  the  total  energy  of  a  relaxed 
discrete  geodesic  between  two  multicomponent  shapes  reads 


{Sl_i)i=l,...,m,  {Sl)i=l,...,m)k=l,...,K] 

^  fl 

=  I  “ ^deform  (‘^fc-l  )i=l,...,m]  +  ^  (^^match  1 ,  +  l^rT^^a  [>5^]  ) 

i=l  \  i=l 

For  sure,  the  different  object  components  or  may  overlap,  but  they  have  to  do  so 
consistently,  that  is,  there  must  exist  a  flow  that  deforms  (Oji)i=i,...,m  into  (CiB)i=i,-  -,m-  In 
fact,  it  is  often  desired  that  the  different  objects  overlap:  Assume  and  to  be  disjoint 
but  have  a  common  boundary.  Obviously,  it  costs  zero  energy  to  pull  both  objects  apart 
rigidly.  Hence,  if  and  0“^  shall  keep  the  common  boundary  along  paths  in  shape  space, 
one  of  the  objects  should  be  replaced  by  the  interior  of  U  O'^  so  that  a  separation  of  both 
components  first  requires  the  costly  generation  of  a  new  boundary.  For  this  reason  we  have 
composed  the  object  in  Fig.  13  of  two  objects,  one  representing  the  whole  bar  and  the  other 
the  black  region.  Another  example  is  given  in  Fig.  14,  where  the  head  and  the  torso  served 
as  one  component  and  the  torso  and  the  legs  as  a  second  one.  Let  us  remark  that  in  case  of 
the  relaxed  model  this  implies  a  different  weighting  of  the  mismatch  penalties  with  respect 
to  different  shape  components. 

Rephrasing  the  above  energy  in  terms  of  level  set  functions  is  straightforward  [42] ,  and 
the  approximations  of  the  different  energy  terms  have  already  been  stated  earlier.  Note 
that  with  m  level  set  functions  and  thus  m  object  components  (9*  we  can  in  fact  distinguish 
n  =  2"*  different  phases  represented  by  objects  0\  i  =  l,...,m,  as  well  as  all  possible 
combinations  of  overlapping.  For  example,  four  phases  (head,  torso,  legs,  background)  have 
been  described  using  two  level  set  functions  in  Fig.  14.  Of  course,  it  is  furthermore  possible 
to  assign  each  phase  different  material  properties.  This  has  been  pursued  in  Fig.  15,  where  a 
geodesic  between  two  frames  from  a  video  of  moving  blood  cells  has  been  computed.  The  top 
row  shows  frames  from  the  real  video  sequence,  where  a  white  blood  cell  squeezes  through  a 
number  of  red  blood  cells.  For  the  computation  of  the  geodesic  (middle  row),  we  employed 
two  level  set  functions  and  assigned  the  white  blood  cell  with  material  parameters  twenty 
times  weaker  than  for  the  red  blood  cells  (material  parameters  of  the  background  are  10^ 
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Figure  14:  Top:  Real  frames  from  a  video  sequence.  Middle:  Discrete  geodesic  between  the 
first  and  the  last  segmented  frames.  Bottom  rows:  Pullback  of  the  last  frame  (top)  and 
pushforward  (bottom)  of  the  first  one  (the  background  has  been  pasted  into  the  pullbacks 
and  pushforwards  so  that  it  is  not  deformed). 


times  weaker).  This  seems  reasonable,  given  that  red  blood  cells  are  comparatively  stiff.  The 
result  is  a  nonlinear  interpolation  between  distant  frames  which  is  in  good  agreement  with 
the  actually  observed  motion.  Note  that  compared  to  Fig.  11  we  need  a  higher  resolution  in 
time  for  Fig.  15  due  to  the  far  stronger  deformations. 

4.3  Application  to  statistics  in  shape  space 

Once  geodesic  distances  between  shapes  are  defined,  one  can  statistically  analyze  ensembles 
of  shapes  and  cluster  them  in  groups  based  on  the  geodesic  distance  as  a  reliable  measure 
for  the  similarity  of  shapes.  A  thorough  investigation  of  such  a  statistical  analysis  is  beyond 
the  scope  of  this  paper.  We  restrict  ourselves  here  to  two  primarily  conceptional  examples  in 
two  and  three  dimensions.  At  first,  we  evalnate  distances  between  different  2D  letters  based 
on  the  discrete  geodesic  path  length.  The  resulting  clustering  is  shown  in  Fig.  16  on  the  left. 
Obviously,  the  Bs  and  As  form  clusters,  and  these  two  clusters  are  closer  to  each  other  than 
the  significantly  distant  M.  Furthermore,  in  Fig.  16  on  the  right  we  study  distances  between 
six  different  foot  level  sets  obtained  from  3D  scans. 

4.4  Complexity  of  shape  space  and  computational  robustness 

In  this  section,  we  intend  to  give  an  impression  of  the  huge  complexity  of  the  Riemannian 
shape  space,  which  is  already  revealed  by  a  small  example.  Fig.  17  shows  a  close-up  of  that 
part  of  shape  space  which  is  spanned  by  the  three  letters  C,  M,  and  U. 

First  note  that  there  may  be  a  large,  possibly  even  infinite  number  of  geodesics  connecting 
these  three  shapes  of  which  we  just  depict  some  interesting  ones.  In  principle,  however,  there 
are  endless  possibilities,  for  example,  to  split  and  remerge  the  letters  in  different  ways  and 
thereby  achieve  a  locally  shortest  path.  Which  path  is  found  by  the  algorithm  depends 
on  the  initialization  of  the  intermediate  shapes  and  deformations  and  therefore  also  on  the 
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Figure  15:  Top:  Frames  from  a  real  video  sequence  of  a  white  blood  cell  among  a  number 
of  red  ones  (courtesy  Robert  A.  Freitas,  Institute  for  Molecular  Manufacturing,  California, 
USA).  Middle:  Computed  discrete  multiphase  geodesic  between  pre- segmented  shapes  in  the 
first  and  the  last  frame.  Bottom:  Pushforward  of  the  initial  (first  four  shapes)  and  pullback 
(last  five  shapes)  of  the  final  frame  according  to  the  geodesic  flow. 
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Figure  16:  Left:  Pairwise  geodesic  distances  between  (also  topologically)  different  letter 
shapes.  Right:  Pairwise  geodesic  distances  between  different  scanned  3D  feet  (data  courtesy 
of  adidas).  The  feet  have  volumes  499.5 cm^,  500.6 cm^,  497.6 cm^,  434.7cm^,  432.0 cm^,  and 
38 1 . 0  cm^ ,  respectively. 
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Figure  17:  Sketch  of  Riemannian  shape  space.  The  box  around  each  geodesic  and  the 
corresponding  path  in  the  sketch  are  coloured  accordingly.  L  denotes  the  geodesic  length 
and  D  the  total  dissipation.  The  circles  represent  the  action  of  on  the  shapes  C,  M,  and 
U  (that  is,  a  rotation  of  the  shapes),  which  induces  no  dissipation  and  has  zero  length. 
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Figure  18:  Discrete  geodesics  between  an  M  and  a  C  of  a  different  font  (cf.  third  row  in 
Fig.  17),  illustrating  robustness  with  respect  to  small  perturbations.  The  geodesic  length 
(and  total  dissipation)  are  0.1220  (0.01518)  for  the  top  row  and  0.1276  (0.01634)  for  the 
bottom  row. 
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Figure  19:  Finer  time  resolution  of  the  geodesics  between  M  and  C  from  Fig.  17  and  18. 
From  top  to  bottom,  the  geodesic  lengths  (and  the  total  dissipation)  are  0.1025  (0.01056), 
0.1201  (0.01465),  and  0.1259  (0.01596). 


position  of  the  end  shapes.  This  position  is  indicated  in  the  sketch  by  the  circle  associated 
to  each  letter,  which  shall  represent  all  possible  rotation  angles. 

It  is  actually  quite  intuitive  that  the  shortest  geodesic  between  the  C  and  the  U  involves 
a  rotation  by  |  (top  geodesic  in  Fig.  17).  Note  that  this  rotation  generates  no  dissipation  so 
that  the  rotated  and  the  upright  C  are  identified  with  each  other  as  being  exactly  the  same 
shape.  Similarly,  the  shortest  geodesics  between  M  and  C  as  well  as  between  M  and  U  are 
such  that  the  inner  two  line  segments  of  the  M  are  bent  outwards  to  yield  a  rotated  C  and 
U  (third  and  last  geodesic).  If  crack  formation  or  closure  is  involved,  the  paths  typically 
exhibit  stronger  dissipation  near  the  cracks  and  are  thus  longer. 

On  the  one  hand  our  computations  are  able  to  reveal  at  least  interesting  parts  of  the 
complex  structure  of  shape  space.  On  the  other  hand,  the  algorithm  proves  to  be  very  robust 
as  discussed  in  the  following.  The  robustness  of  this  type  of  geodesics  becomes  apparent  if 
we  slightly  perturb  the  end  shapes.  Of  course,  we  then  expect  a  similar  geodesic  path  with 
similar  intermediate  shapes,  a  similar  distribution  of  the  dissipation,  and  a  similar  geodesic 
length.  Fig.  18  illustrates  this  continuous  dependence  of  geodesics  on  the  end  shapes,  where 
two  more  versions  of  the  geodesic  between  M  and  C  are  computed. 

As  elaborated  earlier,  our  rigid  body  motion  invariant  time  discretization  is  particularly 
intended  to  allow  good  approximations  to  a  continuous  geodesic  already  for  relatively  coarse 
time  discretizations.  Hence,  as  we  refine  a  discrete  geodesic  by  increasing  the  number  of 
time  steps,  we  expect  that  the  intermediate  shapes  do  not  change  too  strongly  and  that  the 
geodesic  length  has  already  almost  converged.  Indeed,  the  geodesic  length  decreases  by  less 
than  two  percent  when  doubling  K  from  8  to  16  for  the  different  geodesics  between  Ms  and 
Cs  (Fig.  19  and  Table  1). 

We  also  expect  from  a  discrete  geodesic  to  approximately  satisfy  the  axioms  of  a  metric, 
that  is,  that  the  geodesic  distance  between  two  shapes  is  roughly  symmetric  and  satisfies  the 
triangle  inequality.  Let  us  remark  that  in  general  ^deform ^  £deform[4>k^ ,<Sk].  Hence, 
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W  =  4 

iF  =  8 

iF  =  16 

M  O 

M  n 

MO 

0.1068  (0.01145) 
0.1265  (0.01620) 
0.1324  (0.01765) 

0.1040  (0.01084) 
0.1220  (0.01518) 
0.1276  (0.01634) 

0.1025  (0.01056) 
0.1201  (0.01465) 
0.1259  (0.01596) 

Table  1:  Geodesic  distances  (and  dissipation  in  brackets)  for  discrete  geodesics  of  different 
resolutions. 


Figure  20:  Comparison  of  the  discrete  geodesic  from  M  to  C  (copied  from  Fig.  17)  with  the 
one  form  C  to  M.  The  geodesic  lengths  (total  dissipation)  are  0.1040  (0.01084)  and  0.1030 
(0.01064),  respectively.  Note  that  each  shape  from  the  first  row  should  be  compared  to 
the  shape  exactly  below  it,  while  the  color-coding  of  each  shape  in  the  top  row  shows  the 
dissipation  invested  in  the  time  step  to  the  right  whereas  the  color-coding  in  the  bottom  row 
on  each  shape  shows  the  dissipation  in  the  tim  step  to  the  left. 

the  discrete  geodesic  distance  is  by  construction  not  ensured  to  be  strictly  symmetric,  whereas 
the  limit  metric  Q  is  symmetric.  These  two  properties  are  exemplarily  illustrated  in  Fig.  20 
and  21. 


5  Conclusions  and  future  work 

We  have  proposed  a  physical  framework  to  compute  geodesics  in  shape  space  based  on 
measuring  flow-induced  dissipation  in  the  interior  of  shape  contours.  Furthermore,  we  have 
introduced  a  novel  variational  time  discretization  of  these  geodesics.  The  key  ingredients 
are  the  1-1  mapping  property  between  consecutive  time  steps  and  the  rigid  body  motion 
invariance.  The  proposed  formulation  allows  to  weight  the  effect  of  the  local  change  of  length 
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Figure  21:  The  geodesic  lengths  between  the  three  M-  and  C-shapes  satisfy  the  triangle 
inequality  (values  in  brackets  are  total  dissipation). 
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Figure  22:  Limitations  of  the  method:  It  is  advantageous  to  decompose  the  shape  into  small 
independent  pieces  and  to  shuffle  them  around. 


and  volume  separately,  leading  to  significantly  different  geodesic  paths.  Both  physically 
and  with  respect  to  the  shape  description,  geodesic  paths  can  undergo  certain  topological 
transitions.  A  cascadic  multi-scale  relaxation  strategy  renders  the  computation  robust  and 
effective. 

However,  there  are  some  limitations  and  complications  of  our  approach  that  might  be 
tackled  in  future  work.  Since  shapes  are  regarded  as  boundaries  of  volumetric  objects,  we 
need  to  resolve  the  complete  volume  during  computations,  resulting  in  increased  computa¬ 
tion  times  and  thus  restricting  the  method  to  applications  of  relatively  limited  complexity. 
Also,  on  a  classical,  regular  grid  in  space  and  time,  strongly  nonlinear  deformations  such  as 
local  rotations  of  parts  of  an  object  can  be  described  only  at  very  high  resolution,  at  which 
the  problem  relaxation  takes  comparably  long.  Another  issue  concerns  the  perimeter  regu¬ 
larization  without  which  general  minimizing  geodesics  cannot  exist:  It  would  be  optimal  to 
decompose  the  initial  shape  into  a  large  number  of  independent  pieces,  to  shuffle  these  pieces 
around  rigidly,  and  to  recompose  the  end  shape  by  merging  the  pieces  again.  Even  with  reg¬ 
ularization  we  sometimes  observe  a  related  behavior  (see  Fig.  22).  This  problem  is  linked  to 
the  fact  that  there  is  frequently  a  large,  possibly  infinite  number  of  geodesics  between  two 
shapes.  Of  all  possible  geodesics  one  might  select  the  most  intuitive  ones  by  assigning  the 
formation  or  disappearing  of  boundaries  a  further  energy  term.  Finally,  initialization  is  very 
important:  As  we  have  seen  in  Fig.  17,  we  obtain  different  geodesics  depending  on  the  initial 
rotation  angle  of  the  first  shape.  Future  generalizations  of  the  model  might  also  deal  with 
the  incorporation  of  prior  statistical  knowledge.  Furthermore,  we  would  like  to  rigorously 
investigate  the  time-discrete  to  time-continuous  limit  via  the  concept  of  F-convergence. 
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A  Appendix 

A.l  The  Hessian  of  the  deformation  energy 

For  the  deformation  energy  density  W  (X>0) ,  which  we  use  as  an  approximation  to  the  squared 
geodesic  distance  between  two  consecutive  shapes  of  a  discrete  geodesic  path,  there  are  many 
different  possibilities.  As  explained  earlier,  we  only  have  to  ensure  that  W(V(/>)  is  minimal 
and  vanishes  for  G  SO{d),  i.  e.  for  T>(p'^V(p  =  1.  Due  to  W{QT>4>)  =  W{T)(p)  VQ  €  SO{d), 
it  suffices  to  verify  =  0,  where  W^a  denotes  the  derivative  of  W  with  respect  to  its 
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matrix  argument.  Furthermore,  the  second  order  Taylor  expansion  of  W aronnd  0  =  1 
yields  the  desired  local  viscous  dissipation  rate  (1)  as  leading  order  term  for  the  motion  field 
V  =  ^((/)  —  1).  As  an  example — in  fact  used  in  our  computations — we  consider  the  integrand 

^ W)  =  +  ^  (det  V(pf  -  log  det  T>0  -  ^  ^  , 

where  d  G  2,3  is  the  spatial  dimension.  Indeed,  we  obtain  W{1)  =  0  and 

W^aW)  =  i/Dcj)  +  Qdet  T>0  -  cofT>0 , 

which  is  zero  for  Dcj)  =  1,  as  well  as 

W^aaW{T^v,Vv)  =  jiVv  •  +  Q  + 

+  Qdet  1  -  d^)  ^ IVv^]  cofl)  :  Vv 

=  jlDv  :  Vv  +  (/i  +  A)  {tvVvY  —  /i[(trl)u)^  —  Vv^  :  Vv] 

=  2iJ,e[v]  :  e[u]  +  A(tre[u])^  . 


A. 2  Variation  of  the  energy 

Here,  we  give  explicit  formulas  for  the  variation  of  the  different  energy  contributions  in 
directions  of  the  unknown  functions  Uk  {k  =  1, . . . ,  K  —  1)  and  (pk  {k  =  1, ...  ,K),  required 
in  the  numerical  implementation.  Let  ns  denote  by  5^8  {d)  a  variation  of  an  energy  8  with 
respect  to  a  parameter  function  w  in  a  direction  t?.  Using  straightforward  differentiation,  for 
sufficiently  smooth  and  0^  we  obtain 


^0fe^Ltch[0fc,Wfc-i,Wfc](0)  =  2  /  {H^{uk  o  0fc)  -  H^{uk-i))  H'^iuk  o  0fc)(Vufc  o  0^)  •  0da;, 

Jn 

^match 1)  (d)  ~  2  /  (yH^iuk^pk^  Hg(^Uk—i]'&  dx  , 

Jn 

*^«fc^match  '^k  —  ll  ^fc]  (d)  2  /*  (-^£  °  ('Ufc  —  l))  {Uk  O  pk){d  O  pk)  d^  ; 


hu^lLA<t^k,Uk-im  =  I  iil-5)H,{uk-i)  +  S)WAT^^k):Vpdx, 


Su,.,^!Ljcf>k,Uk-im  =  /  il-5)H'{uk-^)^WiVPk)dx, 

Jq 

for  test  functions  t?  and  test  displacements  0,  where  lU^  denotes  the  derivative  of  IF  with 
respect  to  its  matrix  argument.  Employing  the  form  (7)  for  lU  we  obtain 


W,A{A)  =  a,pl!  ^A  +  a2^ll  \oi{A)[l2l  -  cof{A^A)]  +  2T\h)^/hcoi{A) 

V  h 
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for  {Ii,  12,13)  =  (tr(74"'"74),  tr(cof(A"'"74)),  det  . 
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